text:
import .Classification: fit!
using .Classification: GDContext, GradientDescent
using  .NeuralNetwork: NeuralNet, backpropagation, evaluate!

function fit!(net::NeuralNet, context::GDContext, X, y, labels)
    optim = GradientDescent(context)
    res = fit!(net, optim, X, y, labels)
    optim
end

function fit!(net::NeuralNet, optim::GradientDescent, X, y, labels)
    # construct the net
    weave!(net, X, y)
    context = optim.context
    optim.loss = zeros(Float32, context.maxiterations)
    for i in 1:context.maxiterations
        GDs(net, optim, X, y)
        optim.iterations += 1
        evaluate!(net, optim, X, labels)
    end
    postprocess!(optim)
end

function GDs(net::NeuralNet, optimizer::GradientDescent, X::Matrix{<:Real}, y::Matrix{<:Integer})
    N = size(X, 2)
    η = optimizer.context.learningrate
    ∇b = [zeros(Float32, size(layer.b)) for layer in net.layers]
    ∇W = [zeros(Float32, size(layer.W)) for layer in net.layers]
    for sample in 1:size(X, 2)
        b, W = backpropagation(net, X[:, sample],
                               y[:, sample])
        ∇b .+= b
        ∇W .+= W
    end
    for (layer, ∂b, ∂W) in zip(net.layers, ∇b, ∇W)
        layer.b .-= η.*∂b ./N
        layer.W .-= η.*∂W ./N
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_confirmdialogprovider

"""
    dcc_confirmdialogprovider(;kwargs...)
    dcc_confirmdialogprovider(children::Any;kwargs...)
    dcc_confirmdialogprovider(children_maker::Function;kwargs...)

A ConfirmDialogProvider component.
A wrapper component that will display a confirmation dialog
when its child component has been clicked on.

For example:
```
dcc.ConfirmDialogProvider(
    html.Button('click me', id='btn'),
    message='Danger - Are you sure you want to continue.'
    id='confirm')
```
Keyword arguments:
- `children` (Bool | Real | String | Dict | Array; optional): The children to hijack clicks from and display the popup.
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `message` (String; optional): Message to show in the popup.
- `submit_n_clicks` (Real; optional): Number of times the submit was clicked
- `submit_n_clicks_timestamp` (Real; optional): Last time the submit button was clicked.
- `cancel_n_clicks` (Real; optional): Number of times the popup was canceled.
- `cancel_n_clicks_timestamp` (Real; optional): Last time the cancel button was clicked.
- `displayed` (Bool; optional): Is the modal currently displayed.
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
"""
function dcc_confirmdialogprovider(; kwargs...)
        available_props = Symbol[:children, :id, :message, :submit_n_clicks, :submit_n_clicks_timestamp, :cancel_n_clicks, :cancel_n_clicks_timestamp, :displayed, :loading_state]
        wild_props = Symbol[]
        return Component("dcc_confirmdialogprovider", "ConfirmDialogProvider", "dash_core_components", available_props, wild_props; kwargs...)
end

dcc_confirmdialogprovider(children::Any; kwargs...) = dcc_confirmdialogprovider(;kwargs..., children = children)
dcc_confirmdialogprovider(children_maker::Function; kwargs...) = dcc_confirmdialogprovider(children_maker(); kwargs...)

ext: jl

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
"""
    make_unit_cell_FEA_model() -> Dict

This function does, the following in order,

 - Node pairing (from node tags and their coordinate data)
 -

"""
function make_unit_cell_FEA_model(
    ntags::Vector{Int},
    ncoor::Matrix{Float64},
    material_phases_info::Vector{PhaseFiniteElementConnectivity};
    ϵ::Float64 = 1e-06,
    numIP::Int = 4,
    gqOrder::Int = 2,
    verbose::Int = 1,
)::Dict{String, Any}
    #
    # =============================================
    #          MAKING FINITE ELEMENT NODE PAIRS
    # =============================================
    if verbose > 0
        printstyled("FINITE ELEMENT PRE-PROCESSING\n"; color=:yellow, bold=true)
        println(repeat("=", 50))
        println("> Making RVE node pairs,")
    end
    node_pairs = make_rve_node_pairs(ntags, ncoor; small_para = ϵ)
    parent_child_node_tags = [(i.p_nodetag, i.c_nodetag) for i in node_pairs]
    if verbose > 0
        println("> Number of parent-child node pairs ", length(parent_child_node_tags), ",")
    end
    constraining_node_tag = parent_child_node_tags[1][1]  # a_parent_node_tag for constraining DOF
    #
    # =============================================
    #          MAKING FINITE ELEMENT NODE SET
    # =============================================
    if verbose > 0
        println("> Making finite element node sets,")
    end
    all_node_set::Dict{Int, Vector{Float64}} = make_fe_node_set(
        ntags,
        ncoor;
        nodal_dof = ("X", "Y", "Z")
    )
    #
    # =============================================
    #          MAKING FINITE ELEMENT SET
    # =============================================

    element_sets::Vector{FiniteElementSet} = FiniteElementSet[]
    for a_phase_info in material_phases_info
        if verbose > 0
            println("> Making finite element set with ID: '", a_phase_info.tag, "'")
        end
        a_phase_element_sets = make_finite_element_sets(
                a_phase_info.ele_conn,
                all_node_set,
                a_phase_info.material,
                a_phase_info.tag;
                num_ip = numIP,
                gq_order = gqOrder,
            )
        append!(element_sets, a_phase_element_sets)
    end
    if verbose > 0
        println(repeat("=", 50))
    end
    return Dict{String,Any}(
        "all_node_set" => all_node_set,
        "element_sets" => element_sets,
        "pc_node_tags" => parent_child_node_tags,
        "constraing_ntag" => constraining_node_tag,
    )
end


text:
using Reachability, MathematicalSystems, LazySets, TaylorIntegration
using Reachability: solve

# Paremeters
const μ = 3.986e14 * 60^2
const r = 42164.0e3
const r2 = r^2
const mc = 500.0
const n2 = μ / r^3
const n = sqrt(n2)

const two_n = 2*n
const μ_r2 = μ/r2

const K₁ = [-28.8287 0.1005 -1449.9754 0.0046;
            -0.087 -33.2562 0.00462 -1451.5013]
const K₂ = [-288.0288 0.1312 -9614.9898 0.0;
            -0.1312 -288.0 0.0 -9614.9883]

const K₁mc = K₁/mc
const K₂mc = K₂/mc

const zK = zero(K₁)

function switch_controller!(t, x)
    bool_time = constant_term(t) < 120.0
    x_coord = constant_term(constant_term(x[1]))
    bool_approach = -100.0 ≥ x_coord ≥ -1000.0

    if bool_time

        # Approaching mode
        if bool_approach
            return K₁mc * x
        end

        # Rendezvous attempt
        return K₂mc * x
    end

    # Aborting
    return zK * x
end

# dynamics in the 'approaching' mode
@taylorize function space_rendezvous!(t, x, dx)
    x_1 = r + x[1]
    x_12 = x_1^2
    x_22 = x[2]^2
    rc = sqrt(x_12 + x_22)
    rc3 = rc^3
    μ_rc3 = μ / rc3

    uxy = switch_controller!(t, x)

    # x' = vx
    dx[1] = x[3]
    # y' = vy
    dx[2] = x[4]
    # vx' = n²x + 2n*vy + μ/(r^2) - μ/(rc^3)*(r+x) + ux/mc
    dx[3] = (n2*x[1] + two_n*x[4]) + ((μ_r2 - μ_rc3*x_1) + uxy[1])
    # vy' = n²y - 2n*vx - μ/(rc^3)y + uy/mc
    dx[4] = (n2*x[2] - two_n*x[3]) - (μ_rc3*x[2] - uxy[2])
    # t' = 1
    # dx[5] = one(t)

    return dx
end

𝑂 = Options(:T=>200., :plot_vars=>[1, 2], :project_reachset=>true, :mode=>"reach")
X0 = Hyperrectangle([-900., -400., 0., 0.], [25., 25., 0., 0.])
𝑃 = IVP(BlackBoxContinuousSystem(space_rendezvous!, 4), X0)
𝑂jets = Options(:orderT=>10, :orderQ=>2, :abs_tol=>1e-28, :max_steps=>11000)

# solve (long time: 569.110821 seconds (7.69 G allocations: 475.990 GiB, 20.96% gc time)
# @time sol = solve(𝑃, 𝑂, op=TMJets(𝑂jets))

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
module PanZoom

using Gtk
import Gtk.GConstants: GDK_KEY_Left, GDK_KEY_Right, GDK_KEY_Up, GDK_KEY_Down
import Gtk.GConstants.GdkModifierType: SHIFT, CONTROL, MOD1
import Gtk.GConstants.GdkEventMask: KEY_PRESS, SCROLL
import Gtk.GdkEventType: BUTTON_PRESS, DOUBLE_BUTTON_PRESS
import Gtk.GConstants.GdkScrollDirection: UP, DOWN, LEFT, RIGHT
import Base: *, &

using Graphics

using ..GtkUtilities.Link
import ..guidata, ..trigger, ..rubberband_start
import ..Link: AbstractState

const VecLike = Union{AbstractVector,Tuple}

export
    # Types
    Interval,
    # constants
    GDK_KEY_Left,
    GDK_KEY_Right,
    GDK_KEY_Up,
    GDK_KEY_Down,
    SHIFT,
    CONTROL,
    ALT,
    # functions
    interior,
    fullview,
    panzoom,
    panzoom_key,
    panzoom_mouse,
    zoom_reset

const ALT = MOD1

"""
An `Interval` is a `(min,max)` pair. It is the one-dimensional analog
of a `BoundingBox`.
"""
struct Interval{T}
    min::T
    max::T
end

Base.convert(::Type{Interval{T}}, v::Interval{T}) where T = v
Base.convert(::Type{Interval{T}}, v::Interval)  where T = Interval{T}(v.min, v.max)
function Base.convert(::Type{Interval{T}}, v::VecLike) where T
    v1, v2 = first(v), last(v)
    Interval{T}(min(v1,v2), max(v1,v2))
end

Graphics.width(iv::Interval) = iv.max-iv.min
(&)(iv1::Interval, iv2::Interval) = Interval(max(iv1.min, iv2.min),
                                             min(iv1.max, iv2.max))
Graphics.shift(iv::Interval, dx) = deform(iv, dx, dx)
function (*)(s::Real, iv::Interval)
    dx = 0.5*(s - 1)*width(iv)
    deform(iv, -dx, dx)
end
(*)(iv::Interval, s::Real) = s*iv
Graphics.deform(iv::Interval, dmin, dmax) = Interval(iv.min+dmin, iv.max+dmax)

"""
`ivnew = interior(iv, limits)` returns a new version of `iv`, an
`Interval`, which is inside the region allowd by `limits`. One
should prefer "shifting" `iv` over "shrinking" `iv` (if possible, the
width of `iv` should be preserved).

If `limits == nothing`, then `iv` is unconstrained and `ivnew == iv`.

The simplest effectual `limits` object is another `Interval`
representing the full view interval across the chosen axis. If you
need more sophisticated behavior, you can extend this function to work
with custom types of `limits` objects.
"""
interior(iv, ::Nothing) = iv

function interior(iv, limits::Interval)
    imin, imax = iv.min, iv.max
    if imin < limits.min
        imin = limits.min
        imax = imin + width(iv)
    elseif imax > limits.max
        imax = limits.max
        imin = imax - width(iv)
    end
    Interval(imin, imax) & limits
end

"""
`iv = fullview(limits)` returns an `Interval` `iv` that
encompases the full view as permitted by `limits`.

If `limits == nothing`, then `fullview` returns `nothing`.

The simplest effectual `limits` object is another `Interval`
representing the "whole canvas" along the chosen axis. If you need
more sophisticated behavior, you can extend this function to work with
custom types of `limits` objects.
"""
fullview(::Nothing) = nothing

fullview(limits::Interval) = limits

"""
```jl
panzoom(c)
panzoom(c, xviewlimits, yviewlimits)
panzoom(c, xviewlimits, yviewlimits, xview, yview)
panzoom(c2, c1)
```
sets up the Canvas `c` for panning and zooming. The arguments may be
2-tuples, 2-vectors, Intervals, or `nothing`.

`panzoom` creates the `:view[x|y]`, `:view[x|y]limits` properties of
`c`:

- `:xview`, `:yview` are two `AbstractState`s (for horizontal and
    vertical, respectively), each holding an Interval specifying
    the current "view" limits. This might be the entire area, or it
    might be a subregion due to a previous zoom event.

- `:xviewlimits`, `:yviewlimits` encode the maximum allowable viewing
    region; in most cases these will also be `State{Interval}`s, but
    any object that supports `interior` and `fullview` may be used.
    Use `nothing` to indicate unlimited range.

If `c` is the only argument to `panzoom`, then the current user-coordinate
limits of `c` are used.  Note that this invocation works only if the
Canvas has been drawn at least once; if that is not the case, you need
to specify the limits manually.

`panzoom(c2, c1)` sets canvas `c2` to share pan/zoom state with canvas
`c1`.  Panning or zooming in either one will cause the same action in
the other.
"""
panzoom(c, xviewlimits::Interval, yviewlimits::Interval) =
    panzoom(c, State(xviewlimits), State(yviewlimits))

panzoom(c, xviewlimits::VecLike, yviewlimits::VecLike) = panzoom(c, iv(xviewlimits), iv(yviewlimits))

panzoom(c, xviewlimits::Union{VecLike,Nothing}, yviewlimits::Union{VecLike,Nothing}, xview::VecLike, yview::VecLike) = panzoom(c, State(iv(xviewlimits)), State(iv(yviewlimits)), State(iv(xview)), State(iv(yview)))

function panzoom(c, xviewlimits::AbstractState, yviewlimits::AbstractState, xview::AbstractState = similar(xviewlimits), yview::AbstractState = similar(yviewlimits))
    panzoom_disconnect(c)
    if !haskey(guidata, c)
        guidata[c, :xview] = xview
    end
    d = guidata[c]
    d[:xview] = xview
    d[:yview] = yview
    d[:xviewlimits] = xviewlimits
    d[:yviewlimits] = yviewlimits
    link(xview, c)
    link(yview, c)
    nothing
end

const empty_view = State(Interval(0.0, -1.0))

function panzoom(c2, c1)
    panzoom_disconnect(c2)
    d1, d2 = guidata[c1], guidata[c2]
    for s in (:xview, :yview, :xviewlimits, :yviewlimits)
        d2[s] = d1[s]
    end
    link(d2[:xview], c2)
    link(d2[:yview], c2)
    nothing
end

function panzoom(c)
    gc = getgc(c)
    xmin, ymin = device_to_user(gc, 0, 0)
    xmax, ymax = device_to_user(gc, width(c), height(c))
    panzoom(c, (xmin, xmax), (ymin, ymax))
end

function panzoom_disconnect(c)
    for s in (:xview, :yview)
        v = get(guidata, (c, s), empty_view; raw=true)
        if v != empty_view
            disconnect(v, c)
        end
    end
    nothing
end

iv(x) = Interval{Float64}(x...)
iv(x::Interval) = convert(Interval{Float64}, x)
iv(x::State) = iv(get(x))
iv(x::Nothing) = x

pan(iv, frac::Real, limits) = interior(shift(iv, frac*width(iv)), limits)

zoom(iv, s::Real, limits) = interior(s*iv, limits)

"""
`id = panzoom_key(c; kwargs...)` initializes panning- and
zooming-by-keypress for a canvas `c`. `c` is expected to have the four
`guidata` properties described in `panzoom`.

You can configure the keys through keyword arguments. The default
settings are shown below. The first entry is the key, the second a
modifier (like the SHIFT key); `0` means no modifier.

```
    panleft      = (GDK_KEY_Left,0),
    panright     = (GDK_KEY_Right,0),
    panup        = (GDK_KEY_Up,0),
    pandown      = (GDK_KEY_Down,0),
    panleft_big  = (GDK_KEY_Left,SHIFT),
    panright_big = (GDK_KEY_Right,SHIFT),
    panup_big    = (GDK_KEY_Up,SHIFT),
    pandown_big  = (GDK_KEY_Down,SHIFT),
    xpanflip     = false,
    ypanflip     = false
    zoomin       = (GDK_KEY_Up,  CONTROL)
    zoomout      = (GDK_KEY_Down,CONTROL)
```
"Regular" panning motions correspond to 10% of the view region; "big"
panning motions are 100% of the view region, and thus jump by one
whole view area.  The constants are defined in `Gtk.GConstants` and
the modifiers in `Gtk.GConstants.GdkModifierType`.

The returned `id` can be disabled or enabled via
`signal_handler_block` and `signal_handler_unblock`, respectively, or
removed with `signal_handler_disconnect`.

Example:
```
    c = Canvas()
    panzoom(c, (0,1), (0,1))
    id = panzoom_key(c)
```
The `draw` method for `c` should take account of `:viewbb`.
"""
function panzoom_key(c;
                     panleft  = (GDK_KEY_Left,0),
                     panright = (GDK_KEY_Right,0),
                     panup    = (GDK_KEY_Up,0),
                     pandown  = (GDK_KEY_Down,0),
                     panleft_big  = (GDK_KEY_Left,SHIFT),
                     panright_big = (GDK_KEY_Right,SHIFT),
                     panup_big    = (GDK_KEY_Up,SHIFT),
                     pandown_big  = (GDK_KEY_Down,SHIFT),
                     xpanflip     = false,
                     ypanflip     = false,
                     zoomin       = (GDK_KEY_Up,  CONTROL),
                     zoomout      = (GDK_KEY_Down,CONTROL))
    add_events(c, KEY_PRESS)
    set_gtk_property!(c, :can_focus, true)
    set_gtk_property!(c, :has_focus, true)
    signal_connect(key_cb, c, :key_press_event, Cint, (Ptr{Gtk.GdkEventKey},),
                   false, (panleft, panright, panup, pandown, panleft_big,
                           panright_big, panup_big, pandown_big, xpanflip,
                           ypanflip, zoomin, zoomout))
end

@guarded Cint(false) function key_cb(widgetp, eventp, user_data)
    c = convert(GtkCanvas, widgetp)
    event = unsafe_load(eventp)
    (panleft, panright, panup, pandown, panleft_big, panright_big,
     panup_big, pandown_big, xpanflip, ypanflip, zoomin, zoomout) = user_data
    xview = guidata[c, :xview]
    yview = guidata[c, :yview]
    xviewlimits = guidata[c, :xviewlimits]
    yviewlimits = guidata[c, :yviewlimits]
    xsign = xpanflip ? -1 : 1
    ysign = ypanflip ? -1 : 1
    handled = Cint(true)
    ret = Cint(false)
    if keymatch(event, panleft)
        guidata[c, :xview] = pan(xview, -0.1*xsign, xviewlimits)
        ret = handled
    elseif keymatch(event, panright)
        guidata[c, :xview] = pan(xview,  0.1*xsign, xviewlimits)
        ret = handled
    elseif keymatch(event, panup)
        guidata[c, :yview] = pan(yview, -0.1*ysign, yviewlimits)
        ret = handled
    elseif keymatch(event, pandown)
        guidata[c, :yview] = pan(yview,  0.1*ysign, yviewlimits)
        ret = handled
    elseif keymatch(event, panleft_big)
        guidata[c, :xview] = pan(xview, -1*xsign, xviewlimits)
        ret = handled
    elseif keymatch(event, panright_big)
        guidata[c, :xview] = pan(xview,  1*xsign, xviewlimits)
        ret = handled
    elseif keymatch(event, panup_big)
        guidata[c, :yview] = pan(yview, -1*ysign, yviewlimits)
        ret = handled
    elseif keymatch(event, pandown_big)
        guidata[c, :yview] = pan(yview,  1*ysign, yviewlimits)
        ret = handled
    elseif keymatch(event, zoomin)
        xview = zoom(xview, 0.5, xviewlimits)
        yview = zoom(yview, 0.5, yviewlimits)
        setboth(c, xview, yview)
        ret = handled
    elseif keymatch(event, zoomout)
        xview = zoom(xview, 2.0, xviewlimits)
        yview = zoom(yview, 2.0, yviewlimits)
        setboth(c, xview, yview)
        ret = handled
    end
    ret
end

keymatch(event, keydesc) = event.keyval == keydesc[1] && event.state == UInt32(keydesc[2])

"""
`panzoom_mouse(c; kwargs...)` initializes panning-by-mouse-scroll and mouse
control over zooming for a canvas `c`.

zooming or panning (along either x or y) is selected by modifier keys,
which are configurable through keyword arguments.  The default
settings are:

```
    # Panning
    xpan      = SHIFT     # hold down the shift key
    ypan      = 0
    xpanflip  = false
    ypanflip  = false
    # Zooming
    zoom      = CONTROL     # hold down the ctrl-key while scrolling
    focus     = :pointer    # zoom around the position under the mouse pointer
    factor    = 2.0
    initiate  = BUTTON_PRESS # start a rubberband selection for zoom
    reset     = DOUBLE_BUTTON_PRESS    # go back to original limits
```
where 0 means no modifier. SHIFT, CONTROL, BUTTON_PRESS, and
DOUBLE_BUTTON_PRESS are defined in `Gtk.GConstants.GdkModifierType`.

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.
These behaviors are subject to modification by the canvas'
`:viewlimits` data.

An additional keyword is `user_to_data`, for which you may supply
a function
```
    user_to_data_fcn(c, x, y) -> (datax, datay)
```
that converts canvas user-coordinates to "data coordinates" before
setting the values of :xview and :yview.

For important additional information, see `panzoom_key`. To disable mouse
panning and zooming, use
```
    pop!((c.mouse, :scroll))
    pop!((c.mouse, :button1press))
```

Example:
```
    c = Canvas()
    panzoom(c, (0,1), (0,1))
    panzoom_mouse(c)
```
"""
function panzoom_mouse(c;
                       # Panning
                       xpan = SHIFT,
                       ypan  = 0,
                       xpanflip = false,
                       ypanflip  = false,
                       # Zooming
                       zoom = CONTROL,
                       focus::Symbol = :pointer,
                       factor = 2.0,
                       initiate = BUTTON_PRESS,
                       reset = DOUBLE_BUTTON_PRESS,
                       user_to_data = (c,x,y)->(x,y))
    focus == :pointer || focus == :center || error("focus must be :pointer or :center")
    # Scroll events
    scrollfun = (widget, event) -> begin
        s = 0.1*scrollpm(event.direction)
        xscroll = (event.direction == LEFT) || (event.direction == RIGHT)
        if xpan != nothing && (event.state == UInt32(xpan) || xscroll)
            xview = guidata[c, :xview]
            xviewlimits = guidata[c, :xviewlimits]
            guidata[c, :xview] = pan(xview, (xpanflip ? -1 : 1) * s, xviewlimits)
        elseif ypan != nothing && event.state == UInt32(ypan)
            yview = guidata[c, :yview]
            yviewlimits = guidata[c, :yviewlimits]
            guidata[c, :yview] = pan(yview, (ypanflip  ? -1 : 1) * s, yviewlimits)
        elseif zoom != nothing && event.state == UInt32(zoom)
            s = factor
            if event.direction == UP
                s = 1/s
            end
            zoom_focus(widget, s, event; focus=focus, user_to_data=user_to_data)
        end
    end
    # Click events
    clickfun = (widget, event) -> begin
        set_gtk_property!(widget, :is_focus, true)
        if event.event_type == initiate
            rubberband_start(widget, event.x, event.y, (widget, bb) -> zoom_bb(widget, bb, user_to_data))
        elseif event.event_type == reset
            zoom_reset(widget)
        end
    end
    push!((c.mouse, :scroll), scrollfun)
    push!((c.mouse, :button1press), clickfun)
    nothing
end

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

function zoom_focus(c, s, event; focus::Symbol=:pointer, user_to_data=(c,x,y)->(x,y))
    xview = guidata[c, :xview]
    yview = guidata[c, :yview]
    xviewlimits = guidata[c, :xviewlimits]
    yviewlimits = guidata[c, :yviewlimits]
    if focus == :pointer
        ux, uy = device_to_user(getgc(c), event.x, event.y)
        centerx, centery = user_to_data(c, ux, uy)
        w, h = width(xview), width(yview)
        fx, fy = (centerx-xview.min)/w, (centery-yview.min)/h
        wbb, hbb = s*w, s*h
        xview = interior(Interval(centerx-fx*wbb,centerx+(1-fx)*wbb), xviewlimits)
        yview = interior(Interval(centery-fy*hbb,centery+(1-fy)*hbb), yviewlimits)
        setboth(c, xview, yview)
    elseif focus == :center
        xview = zoom(xview, s, xviewlimits)
        yview = zoom(yview, s, yviewlimits)
        setboth(c, xview, yview)
    end
    c
end

function setboth(c, xview, yview)
    getindex(guidata, c, :xview; raw=true).value = xview
    getindex(guidata, c, :yview; raw=true).value = yview
    trigger(c, (:xview, :yview))
end

# We don't take the step of setting new coordinates on the Canvas
# because we need to let the user be in charge of that. (For example,
# in plots you want to zoom in on the data but leave the axes
# visible.) But see set_coordinates below. So we content ourselves with
function zoom_bb(widget, bb::BoundingBox, user_to_data=(c,x,y)->(x,y))
    xmin, ymin = user_to_data(widget, bb.xmin, bb.ymin)
    xmax, ymax = user_to_data(widget, bb.xmax, bb.ymax)
    setboth(widget, (xmin,xmax), (ymin,ymax))
    widget
end

function zoom_reset(widget)
    xvlim, yvlim = guidata[widget, :xviewlimits], guidata[widget, :yviewlimits]
    xvlim != nothing && (getindex(guidata, widget, :xview; raw=true).value = xvlim)
    yvlim != nothing && (getindex(guidata, widget, :yview; raw=true).value = yvlim)
    trigger(widget, (:xview, :yview))
    widget
end

# For completely mysterious reasons, these are borked
# function Graphics.set_coordinates(ctx::GraphicsContext, bb::BoundingBox)
#     set_coordinates(ctx, BoundingBox(0,width(ctx),0,height(ctx)), bb)
# end
# function Graphics.set_coordinates(widget, bb::BoundingBox)
#     set_coordinates(getgc(widget), bb)
# end
function Graphics.set_coordinates(ctx::GraphicsContext, ix::Interval, iy::Interval)
    set_coordinates(ctx, BoundingBox(0,width(ctx),0,height(ctx)), BoundingBox(ix.min, ix.max, iy.min, iy.max))
end

end # module


text:
#=
Blind GridWorld
two observation models are implemented (one is commented out)
Model 1:
The robot is completely blind and observe false all the time
Model 2:
The robot has a noisy measurement of its position, it modeled
by a uniform distribution around the neighbor of the current state (diagonal included)
=#

using StaticArrays
using LinearAlgebra
using Random
using POMDPs
using POMDPModels
using POMDPModelTools
using Compose
using ColorSchemes
using Parameters
using BeliefUpdaters

@with_kw struct BlindGridWorld <: POMDP{GWPos, Symbol, GWPos} #GWPos/Bool
    size::Tuple{Int64, Int64} = (10,10)
    exit::GWPos = GWPos(10,1)
    simple_gw::SimpleGridWorld = SimpleGridWorld(size=size, rewards=Dict(exit=>1.0), terminate_from=Set([exit]))
end

## States

POMDPs.states(pomdp::BlindGridWorld) = states(pomdp.simple_gw)
POMDPs.stateindex(pomdp::BlindGridWorld, s::AbstractVector{Int}) = stateindex(pomdp.simple_gw, s)
POMDPs.initialstate(pomdp::BlindGridWorld) = uniform_belief(pomdp)

## Actions

POMDPs.actions(pomdp::BlindGridWorld) = actions(pomdp.simple_gw)
POMDPs.actionindex(pomdp::BlindGridWorld, a::Symbol) = actionindex(pomdp.simple_gw, a)

## Transition

POMDPs.isterminal(m::BlindGridWorld, s::AbstractVector{Int}) = (s == m.exit || s == GWPos(-1,-1))
POMDPs.transition(pomdp::BlindGridWorld, s::AbstractVector{Int}, a::Symbol) = transition(pomdp.simple_gw, s, a)

## Observation

# Model 1

# POMDPs.observations(pomdp::BlindGridWorld) = (false, true)
# POMDPs.obsindex(pomdp::BlindGridWorld, s::Bool) = Int(s) + 1
# POMDPs.n_observations(pomdp::BlindGridWorld) = 2
# POMDPs.generate_o(pomdp::BlindGridWorld, s, rng) = false

# function POMDPs.observation(pomdp::BlindGridWorld, a::Symbol, sp::AbstractVector{Int})
#     if sp == pomdp.exit
#         return BoolDistribution(0.5)
#     else
#         return BoolDistribution(0.5)
#     end
# end

# Model 2

POMDPs.observations(pomdp::BlindGridWorld) = states(pomdp)
POMDPs.obsindex(pomdp::BlindGridWorld, o) = stateindex(pomdp, o)
function POMDPs.initialobs(pomdp::BlindGridWorld, s)
    ImplicitDistribution(s) do s, rng
        return rand(rng, observation(pomdp, first(actions(pomdp)), s))
    end
end

const NEIGHBORS_DIRECTIONS = Set([GWPos(0,1), GWPos(0,-1), GWPos(-1,0), GWPos(1,0),
                                  GWPos(1,1), GWPos(-1,1), GWPos(1,-1), GWPos(-1,-1)])

function POMDPs.observation(pomdp::BlindGridWorld, a::Symbol, sp::AbstractVector{Int})
    if sp == GWPos(-1,-1)
        return SparseCat(states(pomdp), normalize!(ones(length(states(pomdp))), 1))
    end

    neighbors = MVector{length(NEIGHBORS_DIRECTIONS) + 1, GWPos}(undef)
    neighbors[1] = sp

    # probs = MVector{n_actions(mdp)+1, Float64}()
    probs = @MVector(ones(9))
    for (i, dir) in enumerate(NEIGHBORS_DIRECTIONS)
        neighbors[i+1] = sp + dir

        if !POMDPModels.inbounds(pomdp.simple_gw, neighbors[i+1]) # hit an edge and come back
            probs[i+1] = 0.
        end
    end
    normalize!(probs, 1)

    return SparseCat(neighbors, probs)
end

## Reward
POMDPs.reward(pomdp::BlindGridWorld, s::GWPos) = reward(pomdp.simple_gw, s)
POMDPs.reward(pomdp::BlindGridWorld, s, a, sp) = reward(pomdp.simple_gw, s)
POMDPs.reward(pomdp::BlindGridWorld, s, a) = reward(pomdp.simple_gw, s)
POMDPs.discount(pomdp::BlindGridWorld) = 1.0

## distributions
POMDPs.initialstate(pomdp::BlindGridWorld) = uniform_belief(pomdp)

## helpers

function deterministic_belief(pomdp, s)
    b = zeros(length(states(pomdp)))
    si = stateindex(pomdp, s)
    b[si] = 1.0
    return DiscreteBelief(pomdp, b)
end

function POMDPs.isterminal(pomdp::ProductPOMDP, b::DiscreteBelief)
    belief_support = states(policy.problem)[findall(b.b .> 0)]
    if ProductState(GWPos(8,5),4) ∈ belief_support || ProductState(GWPos(3,7), 3) ∈ belief_support
        return true
    end
    for s in belief_support
        if !isterminal(pomdp, s)
            return false
        end
    end
    return true
end

## Rendering

function POMDPModelTools.render(pomdp::BlindGridWorld, step::Union{NamedTuple, Dict})
    return POMDPModelTools.render(pomdp.simple_gw, step)
end

function POMDPModels.tocolor(r::Float64, minr=0., maxr=1.0)
    frac = (r-minr)/(maxr-minr)
    return get(ColorSchemes.redgreensplit, frac)
end

function POMDPModels.render(mdp::SimpleGridWorld, step::Union{NamedTuple,Dict};
                valuecolor = s -> reward(mdp, s),
                value = nothing,
                action = nothing,
                landmark = nothing,
                minr = 0.,
                maxr = 1.0,
               )

    nx, ny = mdp.size
    cells = []
    vals = []
    acts = []
    landmarks = []
    for x in 1:nx, y in 1:ny
        clr = valuecolor != nothing ? POMDPModels.tocolor(valuecolor(GWPos(x,y)), minr, maxr) : RGB(1.0,1.0,1.0)
        ctx = POMDPModels.cell_ctx((x,y), mdp.size)
        cell = compose(ctx, rectangle(), fill(clr))
        if value != nothing
            val = compose(ctx, text(0.5,0.5, round(value(GWPos(x,y)), digits=3), hcenter, vcenter))
            push!(vals, val)
        end

        if action != nothing
            act = actionarrow(ctx, action(GWPos(x,y)))
            push!(acts, act)
        end

        if landmark != nothing
            if landmark(GWPos(x,y))#haskey(landmark, GWPos(x,y))
                txt = string(LABELLED_STATES[GWPos(x,y)])
                # push!(landmarks, compose(ctx, text(0.5, 0.5, txt, hcenter, vcenter), circle(0.5, 0.5, 0.3), fill("purple")))
                push!(landmarks, compose(ctx, text(0.5, 0.5, txt, hcenter, vcenter)))
                # push!(landmarks, compose(ctx, text(0.5,0.5, landmark[GWPos(x,y)], hcenter, vcenter)))
            end
        end

        push!(cells, cell)
    end
    vals = compose(context(), vals...)
    acts = compose(context(), acts...)
    landmarks = compose(context(), landmarks...)
    grid = compose(context(), linewidth(1Compose.mm), stroke("gray"), cells...)
    outline = compose(context(), linewidth(1Compose.mm), rectangle())

    if haskey(step, :s)
        agent_ctx = POMDPModels.cell_ctx(step[:s], mdp.size)
        agent = compose(agent_ctx, circle(0.5, 0.5, 0.4), fill("orange"))
        if haskey(step, :a)
            agent = compose(agent, actionarrow(agent_ctx, step[:a]))
        end
    else
        agent = nothing
    end

    sz = min(w,h)
    gw_ctx = context((w-sz)/2, (h-sz)/2, sz, sz)
    if value != nothing
        gw_ctx = compose(gw_ctx, vals)
    end
    if action != nothing
        gw_ctx = compose(gw_ctx, acts)
    end
    if landmark != nothing
        gw_ctx = compose(gw_ctx, landmarks)
    end
    return compose(gw_ctx, agent, grid)
end

function actionarrow(ctx, act::Symbol)
    if act == :up
        return compose(ctx, line([(0.5,0.5), (0.5, 0.35)]), stroke(["gray"]), arrow())
    elseif act == :down
        return compose(ctx, line([(0.5,0.5), (0.5, 0.65)]), stroke(["gray"]), arrow())

    elseif act == :right
        return compose(ctx, line([(0.5,0.5), (0.65, 0.5)]), stroke(["gray"]), arrow())

    elseif act == :left
        return compose(ctx, line([(0.5,0.5), (0.35, 0.5)]), stroke(["gray"]), arrow())
    end
end

function highlight_goodstates(ctx, mdp::SimpleGridWorld, states::Vector{GWPos})
    goodstates = []
    for s in states
        x, y = s
        cctx = POMDPModels.cell_ctx((x,y), mdp.size)
        push!(goodstates, compose(cctx, circle(), fill("purple")))
    end
    goodstates = compose(context(), goodstates...)
    compose(goodstates, ctx)
end


text:
# Autogenerated wrapper script for Blosc2_jll for x86_64-apple-darwin
export libblosc2

using Zlib_jll
using Zstd_jll
using Lz4_jll
JLLWrappers.@generate_wrapper_header("Blosc2")
JLLWrappers.@declare_library_product(libblosc2, "@rpath/libblosc2.2.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Zstd_jll, Lz4_jll)
    JLLWrappers.@init_library_product(
        libblosc2,
        "lib/libblosc2.2.0.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module BracketHighlighter

using Tokenize
using Tokenize.Tokens
import Tokenize.Tokens: Token, kind, startpos, endpos, untokenize

using Crayons

import OhMyREPL: add_pass!, PASS_HANDLER

mutable struct BracketHighlighterSettings
    crayon::Crayon
end

const BRACKETMATCHER_SETTINGS =
 BracketHighlighterSettings(Crayon(bold = :true, underline = :true))

function (matcher::BracketHighlighterSettings)(crayons::Vector{Crayon}, tokens::Vector{Token}, cursorpos::Int)
    left_bracket_match, right_bracket_match, matched = bracket_match(tokens, cursorpos)
    !matched && return
    crayons[left_bracket_match] = matcher.crayon
    crayons[right_bracket_match] = matcher.crayon
    return
end

setcrayon!(crayon::Crayon) = BRACKETMATCHER_SETTINGS.crayon = crayon
set_token!(crayon::Crayon) = setcrayon!(crayon)

add_pass!(PASS_HANDLER, "BracketHighlighter", BRACKETMATCHER_SETTINGS, true)

# Takes a string and a cursor index.
# Returns index of left matching bracket, right matching bracket
# and if there was a match at all as a 3 tuple.
const LEFT_DELIMS = [Tokens.LPAREN, Tokens.LSQUARE, Tokens.LBRACE]
const RIGHT_DELIMS = [Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE]
function bracket_match(tokens::Vector{Token}, cursoridx::Int)
    enclosing_token_idx = -1
    char_counter = 0
    # find which token we are currently inside
    for (i, token) in enumerate(tokens)
        char_counter += length(untokenize(token))
        if char_counter >= cursoridx
            enclosing_token_idx = i
            break
        end
    end

    if enclosing_token_idx == -1
        return 0, 0, false
    end
    counts = zeros(Int, length(LEFT_DELIMS))

    left_match, left_idx, right_idx = 0, 0, 0
    # Search to the left
    for idx in enclosing_token_idx:-1:0
        idx == 0 && return 0, 0, false
        token = tokens[idx]
        (i = findfirst(RIGHT_DELIMS, kind(token))) != 0 && (counts[i] -= 1)

        if (i = findfirst(LEFT_DELIMS, kind(token))) != 0
            if counts[i] == 0
                left_match = i
                left_idx = idx
                break
            else
                counts[i] += 1
            end
        end
        idx -= 1
    end

    # Search to the right
    for idx in enclosing_token_idx+1:length(tokens) + 1
        idx == length(tokens) + 1 && return 0, 0, false
        token = tokens[idx]
        LEFT_DELIMS[left_match] == kind(token) && (counts[left_match] -= 1)
        if RIGHT_DELIMS[left_match] == kind(token)
            if counts[left_match] == 0
                return left_idx, idx, true
            else
                counts[left_match] += 1
            end
        end
    end
    return 0, 0 , false
end

end # module


text:
struct MWSL3DIntegrand2{C,O,L}
        test_triangular_element::C
        trial_triangular_element::C
        op::O
        test_local_space::L
        trial_local_space::L
end

function (igd::MWSL3DIntegrand2)(u,v)
        α = igd.op.α
        β = igd.op.β
        γ = igd.op.gamma

        x = neighborhood(igd.test_triangular_element,u)
        y = neighborhood(igd.trial_triangular_element,v)

        r = cartesian(x) - cartesian(y)
        R = norm(r)
        G = exp(-γ*R)/(4π*R)

        f = igd.test_local_space(x)
        g = igd.trial_local_space(y)

        j = jacobian(x) * jacobian(y)

        αjG = α*j*G
        βjG = β*j*G

        G = @SVector[αjG*g[i].value for i in 1:6]
        H = @SVector[βjG*g[i].divergence for i in 1:6]

        @SMatrix[dot(f[i].value,G[j])+f[i].divergence*H[j] for i in 1:6, j in 1:6]
end

function momintegrals!(op::MWSingleLayer3D,
    test_local_space::BDMRefSpace, trial_local_space::BDMRefSpace,
    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

    I, J, K, L = SauterSchwabQuadrature.reorder(
        test_triangular_element.vertices,
        trial_triangular_element.vertices, strat)

    test_triangular_element  = simplex(
        test_triangular_element.vertices[I[1]],
        test_triangular_element.vertices[I[2]],
        test_triangular_element.vertices[I[3]])

    trial_triangular_element = simplex(
        trial_triangular_element.vertices[J[1]],
        trial_triangular_element.vertices[J[2]],
        trial_triangular_element.vertices[J[3]])

    igd = MWSL3DIntegrand2(test_triangular_element, trial_triangular_element,
        op, test_local_space, trial_local_space)
    G = SauterSchwabQuadrature.Sautersauterschwab_parameterized(igd, strat)

    A = levicivita(K) == 1 ? @SVector[1,2] : @SVector[2,1]
    B = levicivita(L) == 1 ? @SVector[1,2] : @SVector[2,1]
    for j ∈ 1:3
        q = L[j]
        for i ∈ 1:3
            p = K[i]
            for n in 1:2
                b = B[n]
                for m in 1:2
                    a = A[m]
                    out[2*(i-1)+m,2*(j-1)+n] += G[2*(p-1)+a,2*(q-1)+b]
                end
            end
        end
    end

    nothing
end

struct MWDL3DIntegrand2{C,O,L}
        test_triangular_element::C
        trial_triangular_element::C
        op::O
        test_local_space::L
        trial_local_space::L
end

function (igd::MWDL3DIntegrand2)(u,v)

        γ = igd.op.gamma

        x = neighborhood(igd.test_triangular_element,u)
        y = neighborhood(igd.trial_triangular_element,v)

        r = cartesian(x) - cartesian(y)
        R = norm(r)
        G = exp(-γ*R)/(4π*R)

        GG = -(γ + 1/R) * G / R * r
        T = @SMatrix [
            0 -GG[3] GG[2]
            GG[3] 0 -GG[1]
            -GG[2] GG[1] 0 ]

        f = igd.test_local_space(x)
        g = igd.trial_local_space(y)

        jx = jacobian(x)
        jy = jacobian(y)
        j = jx*jy

        G1 = j*T*g[1][1]
        G2 = j*T*g[2][1]
        G3 = j*T*g[3][1]

        f1 = f[1][1]
        f2 = f[2][1]
        f3 = f[3][1]

        SMatrix{3,3}(
            dot(f1,G1),
            dot(f2,G1),
            dot(f3,G1),
            dot(f1,G2),
            dot(f2,G2),
            dot(f3,G2),
            dot(f1,G3),
            dot(f2,G3),
            dot(f3,G3),)
end

function momintegrals!(op::MWDoubleLayer3D,
    test_local_space::BDMRefSpace, trial_local_space::BDMRefSpace,
    test_triangular_element, trial_triangular_element, out, strat::SauterSchwabStrategy)

    I, J, K, L = SauterSchwabQuadrature.reorder(
        test_triangular_element.vertices,
        trial_triangular_element.vertices, strat)

    test_triangular_element  = simplex(
        test_triangular_element.vertices[I[1]],
        test_triangular_element.vertices[I[2]],
        test_triangular_element.vertices[I[3]])

    trial_triangular_element = simplex(
        trial_triangular_element.vertices[J[1]],
        trial_triangular_element.vertices[J[2]],
        trial_triangular_element.vertices[J[3]])

    igd = MWDL3DIntegrand2(test_triangular_element, trial_triangular_element,
        op, test_local_space, trial_local_space)
    Q = sauterschwab_parameterized(igd, strat)
    for j ∈ 1:3
        for i ∈ 1:3
            out[i,j] += Q[K[i],L[j]]
        end
    end
    nothing
end


text:
# An implementation of LSMR for the solution of the
# over-determined linear least-squares problem
#
#  minimize ‖Ax - b‖
#
# equivalently, of the normal equations
#
#  A'Ax = A'b.
#
# LSMR is formally equivalent to applying MINRES to the normal equations
# but should be more stable. It is also formally equivalent to CRLS though
# LSMR should be expected to be more stable on ill-conditioned or poorly
# scaled problems.
#
# This implementation follows the original implementation by
# Michael Saunders described in
#
# D. C.-L. Fong and M. A. Saunders, LSMR: An Iterative Algorithm for Sparse
# Least Squares Problems, SIAM Journal on Scientific Computing, 33(5),
# pp. 2950--2971, 2011.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Montreal, QC, May 2015.

export lsmr

"""Solve the regularized linear least-squares problem

  minimize ‖b - Ax‖₂² + λ² ‖x‖₂²

using the LSMR method, where λ ≥ 0 is a regularization parameter.
LSQR is formally equivalent to applying MINRES to the normal equations

  (A'A + λ² I) x = A'b

(and therefore to CRLS) but is more stable.

LSMR produces monotonic residuals ‖r‖₂ and optimality residuals ‖A'r‖₂.
rt is formally equivalent to CRLS, though can be slightly more accurate.

Preconditioners M and N may be provided in the form of linear operators and are
assumed to be symmetric and positive definite. If `sqd` is set to `true`,
we solve the symmetric and quasi-definite system

  [ E   A' ] [ r ]   [ b ]
  [ A  -F  ] [ x ] = [ 0 ],

where E = M⁻¹  and F = N⁻¹.

If `sqd` is set to `false` (the default), we solve the symmetric and
indefinite system

  [ E   A' ] [ r ]   [ b ]
  [ A   0  ] [ x ] = [ 0 ].

In this case, `N` can still be specified and indicates the norm
in which `x` should be measured.
"""
function lsmr(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(A,1)),
              N :: AbstractLinearOperator=opEye(size(A,2)),
              sqd :: Bool=false,
              λ :: Float64=0.0, axtol :: Float64=1.0e-8, btol :: Float64=1.0e-8,
              atol :: Float64=0.0, rtol :: Float64=0.0,
              etol :: Float64=1.0e-8, window :: Int=5,
              itmax :: Int=0, conlim :: Float64=1.0e+8,
              radius :: Float64=0.0, verbose :: Bool=false) where T <: Number

  m, n = size(A)
  size(b, 1) == m || error("Inconsistent problem size")
  verbose && @printf("LSMR: system of %d equations in %d variables\n", m, n)

  # If solving an SQD system, set regularization to 1.
  sqd && (λ = 1.0)
  ctol = conlim > 0.0 ? 1/conlim : 0.0
  x = zeros(T, n)

  # Initialize Golub-Kahan process.
  # β₁ M u₁ = b.
  Mu = copy(b)
  u = M * Mu
  β₁ = sqrt(@kdot(m, u, Mu))
  β₁ == 0.0 && return (x, SimpleStats(true, false, [0.0], [0.0], "x = 0 is a zero-residual solution"))
  β = β₁

  @kscal!(m, 1.0/β₁, u)
  @kscal!(m, 1.0/β₁, Mu)
  Nv = copy(A' * u)
  v = N * Nv
  α = sqrt(@kdot(n, v, Nv))

  ζbar = α * β
  αbar = α
  ρ = 1.0
  ρbar = 1.0
  cbar = 1.0
  sbar = 0.0

  # Initialize variables for estimation of ‖r‖.
  βdd = β
  βd = 0.0
  ρdold = 1.0
  τtildeold = 0.0
  θtilde = 0.0
  ζ = 0.0
  d = 0

  # Initialize variables for estimation of ‖A‖ and cond(A).
  Anorm² = α * α
  maxrbar = 0.0
  minrbar = 1.0e+100;

  # Items for use in stopping rules.
  ctol = conlim > 0.0 ? 1 ./ conlim : 0.0;
  rNorm = β
  rNorms = [rNorm]
  ArNorm = ArNorm0 = α * β
  ArNorms = [ArNorm]

  xENorm² = 0.0
  err_lbnd = 0.0
  err_vec = zeros(T, window)

  verbose && @printf("%5s  %7s  %7s  %7s  %7s  %8s  %8s  %7s\n",
                     "Aprod", "‖r‖", "‖A'r‖", "β", "α", "cos", "sin", "‖A‖²")
  verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e\n",
                     1, β₁, α, β₁, α, 0, 1, Anorm²)

  # A'b = 0 so x = 0 is a minimum least-squares solution
  α == 0.0 && return (x, SimpleStats(true, false, [β₁], [0.0], "x = 0 is a minimum least-squares solution"))
  @kscal!(n, 1.0/α, v)
  @kscal!(n, 1.0/α, Nv)

  h = copy(v)
  hbar = zeros(T, n)

  iter = 0
  itmax == 0 && (itmax = m + n)

  status = "unknown"
  on_boundary = false
  solved = solved_mach = solved_lim = (rNorm <= axtol)
  tired  = iter >= itmax
  ill_cond = ill_cond_mach = ill_cond_lim = false
  zero_resid = zero_resid_mach = zero_resid_lim = false
  fwd_err = false

  while ! (solved || tired || ill_cond)
    iter = iter + 1

    # Generate next Golub-Kahan vectors.
    # 1. βu = Av - αu
    @kscal!(m, -α, Mu,)
    @kaxpy!(m, 1.0, A * v, Mu)
    u = M * Mu
    β = sqrt(@kdot(m, u, Mu))
    if β != 0.0
      @kscal!(m, 1.0/β, u)
      @kscal!(m, 1.0/β, Mu)

      # 2. αv = A'u - βv
      @kscal!(n, -β, Nv)
      @kaxpy!(n, 1.0, A' * u, Nv)
      v = N * Nv
      α = sqrt(@kdot(n, v, Nv))
      if α != 0.0
        @kscal!(n, 1.0/α, v)
        @kscal!(n, 1.0/α, Nv)
      end
    end

    # Continue QR factorization
    (chat, shat, αhat) = sym_givens(αbar, λ)

    ρold = ρ
    (c, s, ρ) = sym_givens(αhat, β)
    θnew = s * α
    αbar = c * α

    ρbarold = ρbar
    ζold = ζ
    θbar = sbar * ρ
    ρtemp = cbar * ρ
    (cbar, sbar, ρbar) = sym_givens(ρtemp, θnew)
    ζ = cbar * ζbar
    ζbar = -sbar * ζbar

    xENorm² = xENorm² + ζ * ζ
    err_vec[mod(iter, window) + 1] = ζ
    iter >= window && (err_lbnd = @knrm2(window, err_vec))

    # Update h, h_har, x.
    hbar = h - (θbar * ρ / (ρold * ρbarold)) * hbar

    # if a trust-region constraint is given, compute step to the boundary
    # the step ϕ/ρ is not necessarily positive
    σ = ζ / (ρ * ρbar)
    if radius > 0.0
      t1, t2 = to_boundary(x, hbar, radius)
      tmax, tmin = max(t1, t2), min(t1, t2)
      on_boundary = σ > tmax || σ < tmin
      σ = σ > 0 ? min(σ, tmax) : max(σ, tmin)
    end

    x = x + σ * hbar
    h = v - (θnew / ρ) * h

    # Estimate ‖r‖.
    βacute =  chat * βdd
    βcheck = -shat * βdd

    βhat =  c * βacute
    βdd  = -s * βacute

    θtildeold = θtilde
    (ctildeold, stildeold, ρtildeold) = sym_givens(ρdold, θbar)
    θtilde = stildeold * ρbar
    ρdold = ctildeold * ρbar
    βd = -stildeold * βd + ctildeold * βhat

    τtildeold = (ζold - θtildeold * τtildeold) / ρtildeold
    τd = (ζ - θtilde * τtildeold) / ρdold
    d = d + βcheck * βcheck
    rNorm = sqrt(d + (βd - τd)^2 + βdd * βdd)
    push!(rNorms, rNorm)

    # Estimate ‖A‖.
    Anorm² += β * β
    Anorm   = sqrt(Anorm²)
    Anorm² += α * α

    # Estimate cond(A).
    maxrbar = max(maxrbar, ρbarold)
    iter > 1 && (minrbar = min(minrbar, ρbarold))
    Acond = max(maxrbar, ρtemp) / min(minrbar, ρtemp)

    # Test for convergence.
    ArNorm = abs(ζbar)
    push!(ArNorms, ArNorm)
    xNorm = @knrm2(n, x)

    test1 = rNorm / β₁
    test2 = ArNorm / (Anorm * rNorm)
    test3 = 1 / Acond
    t1    = test1 / (1.0 + Anorm * xNorm / β₁)
    rNormtol  = btol + axtol * Anorm * xNorm / β₁

    verbose && @printf("%5d  %7.1e  %7.1e  %7.1e  %7.1e  %8.1e  %8.1e  %7.1e\n",
                       1 + 2 * iter, rNorm, ArNorm, β, α, c, s, Anorm²)

    # Stopping conditions that do not depend on user input.
    # This is to guard against tolerances that are unreasonably small.
    ill_cond_mach = (1.0 + test3 <= 1.0)
    solved_mach = (1.0 + test2 <= 1.0)
    zero_resid_mach = (1.0 + t1 <= 1.0)

    # Stopping conditions based on user-provided tolerances.
    tired  = iter >= itmax
    ill_cond_lim = (test3 <= ctol)
    solved_lim = (test2 <= axtol)
    solved_opt = ArNorm <= atol + rtol * ArNorm0
    zero_resid_lim = (test1 <= rNormtol)
    iter >= window && (fwd_err = err_lbnd <= etol * sqrt(xENorm²))

    ill_cond = ill_cond_mach | ill_cond_lim
    solved = solved_mach | solved_lim | solved_opt | zero_resid_mach | zero_resid_lim | fwd_err | on_boundary
  end

  tired         && (status = "maximum number of iterations exceeded")
  ill_cond_mach && (status = "condition number seems too large for this machine")
  ill_cond_lim  && (status = "condition number exceeds tolerance")
  solved        && (status = "found approximate minimum least-squares solution")
  zero_resid    && (status = "found approximate zero-residual solution")
  fwd_err       && (status = "truncated forward error small enough")
  on_boundary   && (status = "on trust-region boundary")

  stats = SimpleStats(solved, !zero_resid, rNorms, ArNorms, status)
  return (x, stats)
end


text:
@testset "add periods" begin
    @test nd + Year(1) == NanoDate(yr + 1, mn, dy, hr, mi, sc, ms, cs, ns)
    @test nd + Month(1) == NanoDate(yr, mn + 1, dy, hr, mi, sc, ms, cs, ns)
    @test nd + Day(1) == NanoDate(yr, mn, dy + 1, hr, mi, sc, ms, cs, ns)
    @test nd + Hour(1) == NanoDate(yr, mn, dy, hr + 1, mi, sc, ms, cs, ns)
    @test nd + Minute(1) == NanoDate(yr, mn, dy, hr, mi + 1, sc, ms, cs, ns)
    @test nd + Second(1) == NanoDate(yr, mn, dy, hr, mi, sc + 1, ms, cs, ns)
    @test nd + Millisecond(1) == NanoDate(yr, mn, dy, hr, mi, sc, ms + 1, cs, ns)
    @test nd + Microsecond(1) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs + 1, ns)
    @test nd + Nanosecond(1) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns + 1)
end

@testset "subtract periods" begin
    @test nd - Year(1) == NanoDate(yr - 1, mn, dy, hr, mi, sc, ms, cs, ns)
    @test nd - Month(1) == NanoDate(yr, mn - 1, dy, hr, mi, sc, ms, cs, ns)
    @test nd - Day(1) == NanoDate(yr, mn, dy - 1, hr, mi, sc, ms, cs, ns)
    @test nd - Hour(1) == NanoDate(yr, mn, dy, hr - 1, mi, sc, ms, cs, ns)
    @test nd - Minute(1) == NanoDate(yr, mn, dy, hr, mi - 1, sc, ms, cs, ns)
    @test nd - Second(1) == NanoDate(yr, mn, dy, hr, mi, sc - 1, ms, cs, ns)
    @test nd - Millisecond(1) == NanoDate(yr, mn, dy, hr, mi, sc, ms - 1, cs, ns)
    @test nd - Microsecond(1) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs - 1, ns)
    @test nd - Nanosecond(1) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns - 1)
end

@testset "trunc to period" begin
    @test trunc(nd, Year) == NanoDate(yr)
    @test trunc(nd, Month) == NanoDate(yr, mn)
    @test trunc(nd, Day) == NanoDate(yr, mn, dy)
    @test trunc(nd, Hour) == NanoDate(yr, mn, dy, hr)
    @test trunc(nd, Minute) == NanoDate(yr, mn, dy, hr, mi)
    @test trunc(nd, Second) == NanoDate(yr, mn, dy, hr, mi, sc)
    @test trunc(nd, Millisecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms)
    @test trunc(nd, Microsecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs)
    @test trunc(nd, Nanosecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns)
end

@testset "floor to period" begin
    @test floor(nd, Year) == NanoDate(yr)
    @test floor(nd, Month) == NanoDate(yr, mn)
    @test floor(nd, Day) == NanoDate(yr, mn, dy)
    @test floor(nd, Hour) == NanoDate(yr, mn, dy, hr)
    @test floor(nd, Minute) == NanoDate(yr, mn, dy, hr, mi)
    @test floor(nd, Second) == NanoDate(yr, mn, dy, hr, mi, sc)
    @test floor(nd, Millisecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms)
    @test floor(nd, Microsecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs)
    @test floor(nd, Nanosecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns)
end

@testset "ceil to period" begin
    @test ceil(nd, Year) == NanoDate(yr + !isone(mn))
    @test ceil(nd, Month) == NanoDate(yr, mn + !isone(dy))
    @test ceil(nd, Day) == NanoDate(yr, mn, dy + !iszero(hr))
    @test ceil(nd, Hour) == NanoDate(yr, mn, dy, hr + !iszero(mi))
    @test ceil(nd, Minute) == NanoDate(yr, mn, dy, hr, mi + !iszero(sc))
    @test ceil(nd, Second) == NanoDate(yr, mn, dy, hr, mi, sc + !iszero(ms))
    @test ceil(nd, Millisecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms + !iszero(cs))
    @test ceil(nd, Microsecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs + !iszero(ns))
    @test ceil(nd, Nanosecond) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns)
end

@testset "round down to period" begin
    @test round(nd, Year, RoundDown) == NanoDate(yr)
    @test round(nd, Month, RoundDown) == NanoDate(yr, mn)
    @test round(nd, Day, RoundDown) == NanoDate(yr, mn, dy)
    @test round(nd, Hour, RoundDown) == NanoDate(yr, mn, dy, hr)
    @test round(nd, Minute, RoundDown) == NanoDate(yr, mn, dy, hr, mi)
    @test round(nd, Second, RoundDown) == NanoDate(yr, mn, dy, hr, mi, sc)
    @test round(nd, Millisecond, RoundDown) == NanoDate(yr, mn, dy, hr, mi, sc, ms)
    @test round(nd, Microsecond, RoundDown) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs)
    @test round(nd, Nanosecond, RoundDown) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns)
end

@testset "round up to period" begin
    @test round(nd, Year, RoundUp) == NanoDate(yr + !isone(mn))
    @test round(nd, Month, RoundUp) == NanoDate(yr, mn + !isone(dy))
    @test round(nd, Day, RoundUp) == NanoDate(yr, mn, dy + !iszero(hr))
    @test round(nd, Hour, RoundUp) == NanoDate(yr, mn, dy, hr + !iszero(mi))
    @test round(nd, Minute, RoundUp) == NanoDate(yr, mn, dy, hr, mi + !iszero(sc))
    @test round(nd, Second, RoundUp) == NanoDate(yr, mn, dy, hr, mi, sc + !iszero(ms))
    @test round(nd, Millisecond, RoundUp) == NanoDate(yr, mn, dy, hr, mi, sc, ms + !iszero(cs))
    @test round(nd, Microsecond, RoundUp) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs + !iszero(ns))
    @test round(nd, Nanosecond, RoundUp) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns)
end

@testset "round nearest (ties up) to period" begin
    @test round(nd, Year, RoundNearestTiesUp) == NanoDate(yr + (mn >= 6))
    @test round(nd, Month, RoundNearestTiesUp) == NanoDate(yr, mn + (dy >= daysinmonth(yr, mn)/2))
    @test round(nd, Day, RoundNearestTiesUp) == NanoDate(yr, mn, dy + (hr >= 12))
    @test round(nd, Hour, RoundNearestTiesUp) == NanoDate(yr, mn, dy, hr + (mn >= 30))
    @test round(nd, Minute, RoundNearestTiesUp) == NanoDate(yr, mn, dy, hr, mi + (sc >= 30))
    @test round(nd, Second, RoundNearestTiesUp) == NanoDate(yr, mn, dy, hr, mi, sc + (ms >= 500))
    @test round(nd, Millisecond, RoundNearestTiesUp) == NanoDate(yr, mn, dy, hr, mi, sc, ms + (cs >= 500))
    @test round(nd, Microsecond, RoundNearestTiesUp) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs + (ns >= 500))
    @test round(nd, Nanosecond, RoundNearestTiesUp) == NanoDate(yr, mn, dy, hr, mi, sc, ms, cs, ns)
end


text:
"""
NodeLoc

A helper struct that points to a node in the tree via its parent such that the child can be easily
swapped out.
If i is 0 the node pointed to is the root node and parent is the node itself.
"""
struct NodeLoc
	parent::RuleNode
	i::Int
end

"""
root_node_loc(root::RuleNode)

Returns a NodeLoc pointing to the root node.
"""
root_node_loc(root::RuleNode) = NodeLoc(root, 0)

"""
get(root::RuleNode, loc::NodeLoc)

Obtain the node pointed to by loc.
"""
function Base.get(root::RuleNode, loc::NodeLoc)
	parent, i = loc.parent, loc.i
	if loc.i > 0
	return parent.children[i]
	else
	return root
	end
end

"""
insert!(loc::NodeLoc, rulenode::RuleNode)

Replaces the subtree pointed to by loc with the given rulenode.
"""
function Base.insert!(root::RuleNode, loc::NodeLoc, rulenode::RuleNode)
	parent, i = loc.parent, loc.i
	if loc.i > 0
		parent.children[i] = rulenode
	else
    		root.ind = rulenode.ind
    		root._val = rulenode._val
    		root.children = rulenode.children
	end
	return root
end


text:
# This file is a part of SimilaritySearch.jl

@inline visited(vstate::BitVector, i)::Bool = @inbounds vstate[i]

@inline function visit!(vstate::BitVector, i)
    @inbounds vstate[i] = true
end

@inline visited(vstate::Vector{UInt8}, i)::Bool = @inbounds vstate[i] == 1

@inline function visit!(vstate::Vector{UInt8}, i, state)
    @inbounds vstate[i] = 1
end

@inline visited(vstate::Set{Int32}, i)::Bool = i ∈ vstate

@inline function visit!(vstate::Set{Int32}, i)
    @inbounds push!(vstate, i)
end

const GlobalVisitedVertices = [BitArray(undef, 1)]  # initialized at __init__ function
# const GlobalVisitedVertices = [Vector{UInt8}(undef, 1)]  # initialized at __init__ function
#const GlobalVisitedVertices = [Set{Int32}()]  # initialized at __init__ function

function __init__visitedvertices()
    for i in 2:Threads.nthreads()
        push!(GlobalVisitedVertices, copy(GlobalVisitedVertices[1]))
    end

    #=for i in eachindex(GlobalVisitedVertices)
        sizehint!(GlobalVisitedVertices[i], 128)
    end=#
end

@inline function _init_vv(v::AbstractVector, n)
    # length(v) < n &&
    resize!(v, n)
    fill!(v, 0)
    v
end

@inline function _init_vv(v::BitVector, n)
    # length(v) < n &&
    resize!(v, n)
    fill!(v.chunks, 0)
    v
end

@inline function _init_vv(v::Set, n)
    empty!(v)
    v
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
using Changepoints
using Base.Test

# write your own tests here
include("test_ChangepointSampler.jl")
include("test_PELT.jl")
include("test_BS.jl")
include("test_R_comparison.jl")


text:
# This file is a part of StruckVMEDevices.jl, licensed under the MIT License (MIT).

import Test
Test.@testset "MemRegisters" begin

include("test_registers.jl")

end # testset


text:
## Return the linear predictor vector
function linpred!(out, p::LinPred, f::Real=1.)
    A_mul_B!(out, p.X, f == 0 ? p.beta0 : broadcast!(muladd, p.scratchbeta, f, p.delbeta, p.beta0))
end
linpred(p::LinPred, f::Real=1.) = linpred!(Vector{eltype(p.X)}(size(p.X, 1)), p, f)

## Install beta0 + f*delbeta as beta0 and zero out delbeta
function installbeta!(p::LinPred, f::Real=1.)
    beta0 = p.beta0
    delbeta = p.delbeta
    @inbounds for i = eachindex(beta0,delbeta)
        beta0[i] += delbeta[i]*f
        delbeta[i] = 0
    end
    p.beta0
end

type DensePredQR{T<:BlasReal} <: DensePred
    X::Matrix{T}                  # model matrix
    beta0::Vector{T}              # base coefficient vector
    delbeta::Vector{T}            # coefficient increment
    scratchbeta::Vector{T}
    qr::QRCompactWY{T}
    function (::Type{DensePredQR{T}}){T}(X::Matrix{T}, beta0::Vector{T})
        n, p = size(X)
        length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
        new{T}(X, beta0, zeros(T,p), zeros(T,p), qrfact(X))
    end
end
(::Type{DensePredQR})(X::Matrix, beta0::Vector) = DensePredQR{eltype(X)}(X, beta0)
convert{T}(::Type{DensePredQR{T}}, X::Matrix{T}) = DensePredQR{T}(X, zeros(T, size(X, 2)))

function delbeta!{T<:BlasReal}(p::DensePredQR{T}, r::Vector{T})
    p.delbeta = p.qr\r
    return p
end

type DensePredChol{T<:BlasReal,C} <: DensePred
    X::Matrix{T}                   # model matrix
    beta0::Vector{T}               # base vector for coefficients
    delbeta::Vector{T}             # coefficient increment
    scratchbeta::Vector{T}
    chol::C
    scratch::Matrix{T}
end
function DensePredChol(X::StridedMatrix)
    F = cholfact!(float(X'X))
    T = eltype(F)
    DensePredChol(AbstractMatrix{T}(X),
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        F,
        similar(X, T))
end

cholpred(X::StridedMatrix) = DensePredChol(X)

cholfactors(c::Cholesky) = c.factors
Base.LinAlg.cholfact!{T<:FP}(p::DensePredChol{T}) = p.chol

Base.LinAlg.cholfact{T<:FP}(p::DensePredQR{T}) = Cholesky{T,typeof(p.X)}(copy(p.qr[:R]), 'U')
function Base.LinAlg.cholfact{T<:FP}(p::DensePredChol{T})
    c = p.chol
    return Cholesky(copy(cholfactors(c)), c.uplo)
end
Base.LinAlg.cholfact!{T<:FP}(p::DensePredQR{T}) = Cholesky{T,typeof(p.X)}(p.qr[:R], 'U')

function delbeta!{T<:BlasReal}(p::DensePredChol{T}, r::Vector{T})
    A_ldiv_B!(p.chol, At_mul_B!(p.delbeta, p.X, r))
    p
end

function delbeta!{T<:BlasReal}(p::DensePredChol{T}, r::Vector{T}, wt::Vector{T})
    scr = scale!(p.scratch, wt, p.X)
    if VERSION < v"0.5.0-dev+4677"
        cholfact!(At_mul_B!(cholfactors(p.chol), scr, p.X), :U)
    else
        cholfact!(Hermitian(At_mul_B!(cholfactors(p.chol), scr, p.X), :U))
    end
    A_ldiv_B!(p.chol, At_mul_B!(p.delbeta, scr, r))
    p
end

type SparsePredChol{T,M<:SparseMatrixCSC,C} <: GLM.LinPred
    X::M                           # model matrix
    Xt::M                          # X'
    beta0::Vector{T}               # base vector for coefficients
    delbeta::Vector{T}             # coefficient increment
    scratchbeta::Vector{T}
    chol::C
    scratch::M
end
function SparsePredChol{T}(X::SparseMatrixCSC{T})
    chol = cholfact(speye(size(X, 2)))
    return SparsePredChol{eltype(X),typeof(X),typeof(chol)}(X,
        X',
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        zeros(T, size(X, 2)),
        chol,
        similar(X))
end

cholpred(X::SparseMatrixCSC) = SparsePredChol(X)

function delbeta!{T}(p::SparsePredChol{T}, r::Vector{T}, wt::Vector{T})
    scr = scale!(p.scratch, wt, p.X)
    XtWX = p.Xt*scr
    c = p.chol = cholfact(Symmetric{eltype(XtWX),typeof(XtWX)}(XtWX, 'L'))
    p.delbeta = c\Ac_mul_B!(p.delbeta, scr, r)
end

Base.cholfact{T}(p::SparsePredChol{T}) = copy(p.chol)
Base.cholfact!{T}(p::SparsePredChol{T}) = p.chol

invchol(x::DensePred) = inv(cholfact!(x))
invchol(x::SparsePredChol) = cholfact!(x) \ eye(size(x.X, 2))
vcov(x::LinPredModel) = scale!(invchol(x.pp), dispersion(x, true))

function cor(x::LinPredModel)
    Σ = vcov(x)
    invstd = similar(Σ, size(Σ, 1))
    for i = eachindex(invstd)
        invstd[i] = 1 / sqrt(Σ[i, i])
    end
    scale!(invstd, scale!(Σ, invstd))
end

stderr(x::LinPredModel) = sqrt.(diag(vcov(x)))

function show(io::IO, obj::LinPredModel)
    println(io, "$(typeof(obj)):\n\nCoefficients:\n", coeftable(obj))
end

ModelFrame(obj::LinPredModel) = obj.fr
ModelMatrix(obj::LinPredModel) = obj.pp.X
model_response(obj::LinPredModel) = obj.rr.y

fitted(m::LinPredModel) = m.rr.mu
formula(obj::LinPredModel) = ModelFrame(obj).formula
residuals(obj::LinPredModel) = residuals(obj.rr)
predict(mm::LinPredModel) = fitted(mm)
predict(mm::LinPredModel, interval_type::Symbol, level::Real = 0.95) = predict(mm, mm.pp.X, interval_type, level)

# an internal function called by predict for both LinearModel and AbstractGLM
function predict_confint(mm::LinPredModel, newx::AbstractMatrix, interval_type::Symbol, level)
    retmean = newx * coef(mm)
    R = cholfact!(mm.pp)[:U] #get the R matrix from the QR factorization
    residvar = (ones(size(newx,2),1) * deviance(mm)/dof_residual(mm))
    retvariance = (newx/R).^2 * residvar

    interval = quantile(TDist(dof_residual(mm)), (1 - level)/2) * sqrt.(retvariance)
    hcat(retmean, retmean .+ interval, retmean .- interval)
end

"""
    nobs(obj::LinearModel)
    nobs(obj::GLM)

For linear and generalized linear models, returns the number of rows, or,
when prior weights are specified, the sum of weights.
"""
function nobs(obj::LinPredModel)
    if isempty(obj.rr.wts)
        oftype(sum(one(eltype(obj.rr.wts))), length(obj.rr.y))
    else
        sum(obj.rr.wts)
    end
end

coef(x::LinPred) = x.beta0
coef(obj::LinPredModel) = coef(obj.pp)

dof_residual(obj::LinPredModel) = nobs(obj) - length(coef(obj))


text:
neid_data_path = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/"
proc_version_path = "DRPmaster/solar"
path_start = joinpath(neid_data_path,proc_version_path)
output_path_base = "/gpfs/group/ebf11/default/ebf11/neid_solar/data/outputs_1.0.0"
proj_dir = "/gpfs/group/ebf11/default/ebf11/neid_solar/code/NeidSolarScripts.jl"
script = joinpath(proj_dir,"scripts","make_manifest_solar_1.0.0.jl")

sleep_interval = 5

using CSV, DataFrames

dirs_to_process_file = joinpath(output_path_base,"dirs_to_process.csv")
if isfile(dirs_to_process_file)
   println("# Reusing ", dirs_to_process_file)
   df = CSV.read(dirs_to_process_file, DataFrame)
else
   df = DataFrame(:input_dir=>String[],:output_dir=>String[])

   for (root, dirs, files) in walkdir(path_start)
      if !occursin("level1",root)  continue end
      for dir in dirs
              if !occursin(dir,root) continue end
              println(dir) # path to directories
              #println(joinpath(root, dir)) # path to directories
              output_dir = joinpath(output_path_base,dir)
              if !(isdir(output_dir)||islink(output_dir))
                 mkdir(output_dir)
              end
              #continuum_dir = joinpath(output_dir,"continuum")
              #if !(isdir(continuum_dir)||islink(continuum_dir))
              #   mkdir(continuum_dir)
              #end
              push!(df,(input_dir=joinpath(root,dir), output_dir=output_dir))
      end
      #=
      println("Files in $root")
      for file in files
          println(joinpath(root, file)) # path to files
      end
      =#
      # break
   end
   CSV.write(joinpath(output_path_base,"dirs_to_process.csv"),df)
end

last_job_id = 0
function gen_job_name(; prefix::String = "auto")
  global last_job_id
  last_job_id += 1
  return prefix * string(last_job_id)
end

function make_pbs_scr(;proj_dir::String, script::String, input_dir::String, output_dir::String, job_name::String = gen_job_name() )
   pbs_scr = """
#!/bin/bash
#PBS -N $job_name
#PBS -l nodes=1:ppn=1
#PBS -l pmem=2000mb
#PBS -l walltime=4:00:00
###PBS -A ebf11_c_g_vc_default
###PBS -q hprc
#PBS -A cyberlamp
#PBS -l feature=rhel7
#PBS -j oe
#PBS -M ebf11@psu.edu

# Get started
echo Job started on `hostname` at `date`

# Go to the correct place
#cd \$PBS_O_WORKDIR
cd $proj_dir

# Run the job itself
~/julia1.6 --project=$proj_dir -e 'target_subdir=\"$input_dir\"; output_dir=\"$output_dir\";  include(\"$script\")'

# Finish up
echo Job Ended at `date`
"""
  return pbs_scr
end

open("submit_make_manifests.sh","w") do f_submit

for row in eachrow(df)
   # #=
    if isfile(joinpath(row.output_dir, "manifest.csv"))
       println("# Skipping ", row.output_dir)
       continue
    end
   # =#
    m = match(r"(\d+)$", row.input_dir)
    dir = m.captures[1]

    pbs_scr = make_pbs_scr(proj_dir=proj_dir, script=script, input_dir=row.input_dir, output_dir=row.output_dir, job_name = "manif_"*dir)

    #println("Created script: ")
    #println(pbs_scr)
    #println("Echoing script: ")
    open("make_manifest_$dir.pbs","w") do f_pbs
       print(f_pbs, pbs_scr)
    end
    println(f_submit, "qsub make_manifest_$dir.pbs")
    println(f_submit, "sleep $sleep_interval")
    #output = readchomp(pipeline(`echo $pbs_scr`, `qsub`))
    #println("output = ")
    #println(output)
end

end # submit_make_manifests.sh


text:
# # Temperature dependencies

## load packages
using Land.Photosynthesis
using PlotPlants
FT = Float32;
#------------------------------------------------------------------------------

# ## Arrhenius correction
# ### Without deactivation term
## KcTDCLM is a ArrheniusTD type struct, use it as an example here
_td = KcTDCLM(FT);
_ts = collect(FT, 273:1:323);
_ks = temperature_correction.([_td], _ts);

_fig,_axes = create_canvas("Arrhenius correction without deactivation");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _ks, "k-");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Relative to 25 °C");
_fig
#------------------------------------------------------------------------------

# ### With deactivation term
## VcmaxTDCLM is a ArrheniusPeakTD type struct, use it as an example here
_td = VcmaxTDCLM(FT);
_ts = collect(FT, 273:1:323);
_ks = temperature_correction.([_td], _ts);

_fig,_axes = create_canvas("Arrhenius correction with deactivation");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _ks, "k-");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Relative to 25 °C");
_fig
#------------------------------------------------------------------------------

# ## Q10 correction
_td = Q10TD{FT}(1.0, 273.15, 1.7);
_ts = collect(FT, 273:1:323);
_ks = temperature_correction.([_td], _ts);

_fig,_axes = create_canvas("Q10 correction");
_ax1 = _axes[1];
_ax1.plot(_ts .- 273.15, _ks, "k-");
_ax1.set_xlabel("Leaf temperature (°C)");
_ax1.set_ylabel("Relative to 25 °C");
_fig
#------------------------------------------------------------------------------

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module Tst

using Speckles

function run()
    # plotsPath = mkpath("plots/20210531_natoms")
    # dataPath = mkpath("data/20210531_natoms")
    plotsPath = mkpath("plots/20210621_roberto1")
    dataPath = mkpath("data/20210621_roberto1")
    # specify all parameters
    # νHα1 = [456808] #GHz
    νHα2 = [456810,456813] #GHz
    # νHα = [456812, 456808,456811, 456802] #GHz

    paramDict = Dict(
                    :n    => [100],#,20,40,80,160], # number of atoms
                    :νm   => [νHα2], # line frequencies in GHz
                    :Em   => ["ones"], # relative line magnitudes
                    :σ    => [20.0], # Doppler broadening in GHz
                    :fγ   => [2.0e6,"shot10%","shot50%",10.0,1.0,0.16], # mean photon count rate in GHz
                    :deadtime   => [0.0], # detector deadtime in nanoseconds
                    :resolution => [0.010],#,0.10], # detector resolution in nanoseconds
                    :jitter     => [0.015], # detector timing jitter in nanoseconds
                    :efficiency => [0.9], # detector efficiency
                    :darkcounts => [1.0e-8], # detector dark count rate in GHz
                    :duration   => [100.0], # duration of each correlation measurement in nanoseconds
                    :repeat     => [100], # number of times to repeat correlation measurement
                    :reinstance => [true], # control whether or not frequencies and phases should be reinstanced between measurements
                    :timeint    => ["halfwindow"] # time over which to average correlations in nanoseconds
                    # :directory  => [] # defaults to main package directory
                    )

    # split into vector of dictionaries: one for each run
    iterParams = paramVector(paramDict)

    # define beamsplitter
    bs = Beamsplitter(0.5,0.5)

    # iterate over parameter dictionaries
    for params in iterParams
        if params[:Em] == "ones" # sets all line magnitudes equal to one if true
            params[:Em] = ones(length(params[:νm]))
        end
        if params[:fγ] == "shot1%"
            params[:fγ] = 2*1e4/params[:resolution] # multiply by 2 so error in each beam is ~1%
        elseif params[:fγ] == "shot10%"
            params[:fγ] = 2*1e2/params[:resolution] # multiply by 2 so error in each beam is ~10%
        elseif params[:fγ] == "shot50%"
            params[:fγ] = 2*4/params[:resolution] # multiply by 2 so error in each beam is ~50%
        end
        if params[:timeint] == "halfwindow"
           params[:timeint] = params[:duration]/2
        end

        # create naming prefix
        prefix = makeName(params)

        # create LightSource and Detector objects
        source = LightSource(
            params[:n],
            params[:Em],
            params[:νm],
            params[:σ],
            params[:fγ]
        )
        detect = Detector(
            params[:deadtime],
            params[:resolution],
            params[:jitter],
            params[:efficiency],
            params[:darkcounts]
        )

        # calculate the average photon counts for each time bin
        γint = γIntensity(
            nbar(params[:duration],source)*bs.t, # apply beamsplitter
            params[:duration],
            detect.resolution,
            source
            )

        # calculate the length in indices of the time integration window
        indexint = length(γint)÷convert(Int,ceil(params[:duration]/params[:timeint]))

        # generate the correlation time offsets
        τ = params[:resolution]*collect(0:(indexint-1))

        # create a 5 e-fold low cut in τ to get rid of non-periodic part of g2
        τCutLow = sqrt(5)/params[:σ] # low τ cut in nanoseconds

        # calculate index of τ cut
        iτCutLow = convert(Int,ceil(τCutLow/params[:resolution]))

        # get the frequencies for a Fourier transform over τ with cuts
        freqs = fftFreq(params[:resolution],τ,(iτCutLow,lastindex(τ)))

        # initialize time and frequency DataFrames
        timeData = DataFrame(:time=>τ)
        freqData = DataFrame(:freq=>freqs[freqs .>= 0])

        # iterate over desired number of repeats
        for i=1:params[:repeat]
            # read out each beam
            readout1 = denseReadout(γint,detect)
            readout2 = denseReadout(γint,detect)

            # plot intensity for the first instance
            if i == 1
                times = params[:resolution]*collect(0:(length(γint)-1))
                # γIntensityPlot(times,γint.γvec,plotsDir(prefix,plotsPath))
                γCountPlot(times .+ params[:resolution]/2,readout1,γint.γvec,plotsDir(prefix,plotsPath))
            end
            # calculate correlation
            corr12 = map(offset->ncorrelate(readout1,readout2,offset,indexint),0:(indexint-1))
            # store instance in dataframe
            iname = "corr$i"
            timeData[!,iname] = corr12

            # take the Fourier transform of the correlation function
            ftcorr12 = meanFFT(corr12,(iτCutLow,lastindex(corr12)))
            ftcorr12pos, = fftPositiveFreq(ftcorr12,freqs)

            freqData[!,iname] = abs.(ftcorr12pos)

            # generate new intensity instance if desired
            if params[:reinstance]
                γint = γIntensity(
                    nbar(params[:duration],source)*bs.t, # multiply by 0.5 for two beams
                    params[:duration],
                    detect.resolution,
                    source
                    )
            end
        end

        # create sum column in frequency and time DataFrames
        timeData[!,:sum] = sum(eachcol(timeData)[2:end])
        freqData[!,:sum] = sum(eachcol(freqData)[2:end])

        # make names for data files
        timeDataName = dataDir(string(prefix,"g2-vs-tau.csv"),dataPath)
        freqDataName = dataDir(string(prefix,"ftg2-vs-freq.csv"),dataPath)

        # save data
        CSV.write(timeDataName,timeData)
        CSV.write(freqDataName,freqData)

        # create and save plots
        γCorrTimePlot(timeData, params, plotsDir(prefix,plotsPath))
        γCorrFreqPlot(freqData, params, plotsDir(prefix,plotsPath))
    end
    return nothing
end #run
export run
end #Tst
import .Tst

Tst.run()


text:
#########################################################################
#
#   ExNode type definition and related functions
#
#########################################################################

#####  ExNode type  ######

type ExNode{T}
  main                    # main payload
  parents::Vector{Any}    # parent nodes
  precedence::Vector{Any} # nodes that should be evaluated before (but are not parents)
  val                     # value
  alloc::Bool             # Allocation ? Forbids fusions

  ExNode(main)                           = new(   main,   Any[], Any[], NaN, false)
  ExNode(main,parents)                   = new(   main, parents, Any[], NaN, false)
  ExNode(main,parents, prec, val, alloc) = new(   main, parents,  prec, val, alloc)
end

copy{T}(x::ExNode{T}) = ExNode{T}( x.main, # copy(x.main),
                                  copy(x.parents),
                                  copy(x.precedence),
                                  # copy(x.val),
                                  x.val,
                                  x.alloc)

copy(x::ExNode{:for}) = ExNode{:for}(Any[ x.main[1], copy(x.main[2]) ],    # make a copy of subgraph
                              copy(x.parents),
                              copy(x.precedence),
                              # copy(x.val),
                              x.val,
                              x.alloc)

typealias NConst     ExNode{:constant}  # for constant
typealias NExt       ExNode{:external}  # external var
typealias NCall      ExNode{:call}      # function call
typealias NComp      ExNode{:comp}      # comparison operator
typealias NRef       ExNode{:ref}       # getindex
typealias NDot       ExNode{:dot}       # getfield
typealias NSRef      ExNode{:subref}    # setindex
typealias NSDot      ExNode{:subdot}    # setfield
typealias NFor       ExNode{:for}       # for loop
typealias NIn        ExNode{:within}    # reference to var set in a loop

subtype{T}(n::ExNode{T}) = T

function show(io::IO, res::ExNode)
  pl = join( map(x->isa(x,NFor) ? "subgraph" : repr(x.main), res.parents) , " / ")
  print(io, "[$(subtype(res))] $(repr(res.main)) ($(repr(res.val)))")
  length(pl) > 0 && print(io, ", from = $pl")
end


text:
using Ipopt
using Compat
using Compat.Test

@testset "C API" begin
    # First of all, test that hs071 example works
    include("hs071_test.jl")

    # Test setting some options
    # String option
    println("\nString option")
    addOption(prob, "hessian_approximation", "exact")
    @test_throws ErrorException addOption(prob, "hessian_approximation", "badoption")
    println("\nInt option")
    # Int option
    addOption(prob, "file_print_level", 3) == nothing
    @test_throws ErrorException addOption(prob, "file_print_level", -1)
    # Double option
    println("\nFloat option")
    addOption(prob, "derivative_test_tol", 0.5)
    @test_throws ErrorException addOption(prob, "derivative_test_tol", -1.0)

    # Test opening an output file
    openOutputFile(prob, "blah.txt", 5)
    Ipopt.freeProblem(prob) # Needed before the `rm` on Windows.
    # unlink the output file
    rm("blah.txt")

    # Test that the ipopt binary works
    # See https://github.com/JuliaOpt/Ipopt.jl/issues/119 for discussion of the
    # known failure on Windows and Julia 0.7.
    if !(Compat.Sys.iswindows() && VERSION >= v"0.7-")
        @test success(`$(Ipopt.amplexe) -v`)
    end
end

@testset "MathProgBase" begin
    include("MPBWrapper.jl")
end

@testset "MathOptInterface" begin
    include("MOIWrapper.jl")
end


text:
# Ref https://arxiv.org/pdf/1901.10739.pdf
# This term does not contain the kinetic energy, which must be added separately
# /!\ They have no 1/2 factor in front of the kinetic energy,
#     so for consistency the added kinetic energy must have a scaling of 2
# Energy = <u, ((-i∇ + βA)^2 + V) u>
# where ∇∧A = 2π ρ, ∇⋅A = 0 => A = x^⟂/|x|² ∗ ρ
# H = (-i∇ + βA)² + V - 2β x^⟂/|x|² ∗ (βAρ + J)
#   = -Δ + 2β (-i∇)⋅A + β²|A|^2 - 2β x^⟂/|x|² ∗ (βAρ + J)
# where only the first three terms "count for the energy", and where
# J = 1/(2i) (u* ∇u - u ∇u*)

struct Anyonic
    β
end
function (A::Anyonic)(basis)
    @assert length(basis.kpoints) == 1
    @assert basis.kpoints[1].coordinate == [0, 0, 0]
    @assert basis.model.dim == 2
    # only square lattices allowed
    # (because I can't be bothered to think about the right formulas otherwise,
    # although they might already be correct)
    @assert basis.model.lattice[2, 1] == basis.model.lattice[1, 2] == 0
    @assert basis.model.lattice[1, 1] == basis.model.lattice[2, 2]

    TermAnyonic(basis, eltype(basis)(A.β))
end

struct TermAnyonic{T <: Real} <: Term
    basis::PlaneWaveBasis{T}
    β::T
end

function ene_ops(term::TermAnyonic, ψ, occ; ρ, kwargs...)
    basis = term.basis
    T = eltype(basis)
    β = term.β
    @assert ψ !== nothing # the hamiltonian depends explicitly on ψ

    # Compute A in Fourier domain
    # curl A = 2π ρ, ∇⋅A = 0
    # i [G1;G2] ∧ A(G1,G2) = 2π ρ(G1,G2), [G1;G2] ⋅ A(G1,G2) = 0
    # => A(G1, G2) = -2π i ρ(G1, G2) * [-G2;G1;0]/(G1^2 + G2^2)
    A1 = zeros(complex(T), basis.fft_size)
    A2 = zeros(complex(T), basis.fft_size)
    for (iG, G) in enumerate(G_vectors_cart(basis))
        G2 = sum(abs2, G)
        if G2 != 0
            A1[iG] = +2T(π) * G[2] / G2 * ρ.fourier[iG] * im
            A2[iG] = -2T(π) * G[1] / G2 * ρ.fourier[iG] * im
        end
    end
    Areal = [from_fourier(basis, A1).real,
             from_fourier(basis, A2).real,
             zeros(T, basis.fft_size)]

    # 2β (-i∇)⋅A + β^2 |A|^2
    ops_energy = [MagneticFieldOperator(basis, basis.kpoints[1],
                                        2β .* Areal),
                  RealSpaceMultiplication(basis, basis.kpoints[1],
                                          β^2 .* (abs2.(Areal[1]) .+ abs2.(Areal[2])))]

    # Now compute effective local potential - 2β x^⟂/|x|² ∗ (βAρ + J)
    J = compute_current(basis, ψ, occ)
    eff_current = [J[α].real .+
                   β .* ρ.real .* Areal[α] for α = 1:2]
    eff_current_fourier = [from_real(basis, eff_current[α]).fourier for α = 1:2]
    # eff_pot = - 2β x^⟂/|x|² ∗ eff_current
    # => ∇∧eff_pot = -4πβ eff_current
    # => eff_pot(G1, G2) = 4πβ i eff_current(G1, G2) * [-G2;G1;0]/(G1^2 + G2^2)
    eff_pot_fourier = zeros(complex(T), basis.fft_size)
    for (iG, Gred) in enumerate(G_vectors(basis))
        G = basis.model.recip_lattice * Gred
        G2 = sum(abs2, G)
        if G2 != 0
            eff_pot_fourier[iG] += -4T(π)*β * im * G[2] / G2 * eff_current_fourier[1][iG]
            eff_pot_fourier[iG] += +4T(π)*β * im * G[1] / G2 * eff_current_fourier[2][iG]
        end
    end
    eff_pot_real = from_fourier(basis, eff_pot_fourier).real
    ops_ham = [ops_energy..., RealSpaceMultiplication(basis, basis.kpoints[1], eff_pot_real)]

    E = zero(T)
    for iband = 1:size(ψ[1], 2)
        ψnk = @views ψ[1][:, iband]
        # TODO optimize this
        for op in ops_energy
            E += occ[1][iband] * real(dot(ψnk, op * ψnk))
        end
    end

    (E=E, ops=[ops_ham])
end


text:
export GADM, NUTS, makeregions, makeregions_nuts, makeoffshoreregions, saveregions, loadregions,
        saveregions_global, subregions, all_gadm_subregions

abstract type RegionType end

struct GADM{T} <: RegionType
    parentregions::Vector{T}
    subregionnames::NTuple{N,T} where N
end
GADM(regionnames::T...) where T = GADM(T[], regionnames)
GADM(parentregions::Vector{T}, subregionnames::T...) where T = GADM(parentregions, subregionnames)

struct NUTS{T} <: RegionType
    subregionnames::NTuple{N,T} where N
end
NUTS(regionnames::T...) where T = NUTS(regionnames)

const NOREGION = typemax(Int16)

function saveregions(regionname, regiondefinitionarray; autocrop=true, bbox=[-90 -180; 90 180])
    land = JLD.load(in_datafolder("landcover.jld"), "landcover")
    if !all(bbox .== [-90 -180; 90 180])
        autocrop = false         # ignore supplied autocrop option if user changed bbox
    end
    saveregions(regionname, regiondefinitionarray, land, autocrop, bbox)
end

function saveregions(regionname, regiondefinitionarray, landcover, autocrop, bbox)
    regions, regiontype = makeregions(regiondefinitionarray; allowmixed=(regionname=="Europe_background"))
    if autocrop
        # get indexes of the bounding box containing onshore region data with 6% of padding
        lonrange, latrange = getbboxranges(regions)
        padding = round(Int, maximum(size(regions[lonrange,latrange])) * 0.06)
        lonrange, latrange = getbboxranges(regions, padding)
    else
        latrange, lonrange = bbox2ranges(roundbbox(bbox,100), 100)          # TO DO: remove hardcoded raster density
    end
    landcover = landcover[lonrange, latrange]
    regions = regions[lonrange, latrange]

    if regionname != "Global_GADM0" && regionname != "Europe_background"
        if regiontype == :NUTS
            println("\nNUTS region definitions detected (using Europe_background region file)...")
            europeregions = loadregions("Europe_background")[1][lonrange, latrange]
            regions[(regions.==0) .& (europeregions.>0)] .= NOREGION
        elseif regiontype == :GADM
            println("\nGADM region definitions detected (using Global_GADM0 region file)...")
            globalregions = loadregions("Global_GADM0")[1][lonrange, latrange]
            regions[(regions.==0) .& (globalregions.>0)] .= NOREGION
        end
    end

    # Find the closest region pixel for all non-region pixels (land and ocean)
    println("\nAllocate non-region pixels to the nearest region (for offshore wind)...")
    territory = regions[feature_transform(regions.>0)]

    # Allocate ocean and lake pixels to the region with the closest land region.
    # Even VERY far offshore pixels will be allocated to whatever region is nearest, but
    # those areas still won't be available for offshore wind power because of the
    # requirement to be close enough to the electricity grid (or rather the grid proxy).
    offshoreregions = territory .* (landcover .== 0)

    if regionname != "Global_GADM0" && regionname != "Europe_background"
        # Allocate land pixels with region==0 to the closest land region.
        # This ensures that the regions dataset is pixel-compatible with the landcover dataset.
        regions = territory .* (landcover .> 0)
    end

    println("\nSaving regions and offshoreregions...")
    regionlist = Symbol.(regiondefinitionarray[:,1])

    JLD.save(in_datafolder("regions_$regionname.jld"), "regions", regions, "offshoreregions", offshoreregions,
                "regionlist", regionlist, "lonrange", lonrange, "latrange", latrange, compress=true)
end

function saveregions_global(; args...)
    println("\nCreating a global GADM region file to identify countries and land areas later...")
    g = readdlm(in_datafolder("gadmfields.csv"), ',', skipstart=1)
    gadm0 = unique(string.(g[:,2]))
    regiondefinitionarray = [gadm0 GADM.(gadm0)]
    saveregions("Global_GADM0", regiondefinitionarray; args..., autocrop=false)
    println("Global GADM region file saved.")

    println("\nCreating a 'background' NUTS region file to identify non-European land areas later...")
    regiondefinitionarray = [NUTS_Europe; non_NUTS_Europe]
    saveregions("Europe_background", regiondefinitionarray; args..., autocrop=false)
    println("\nEurope_background region file saved.")

    println("\nCreating a region file for the 44 countries with synthetic demand training data...")
    regiondefinitionarray = [syntheticdemandregions GADM.(syntheticdemandregions)]
    saveregions("SyntheticDemandRegions", regiondefinitionarray; args...)
    println("\nSyntheticDemandRegions file saved.")
end

function loadregions(regionname)
    jldopen(in_datafolder("regions_$regionname.jld"), "r") do file
        return read(file, "regions"), read(file, "offshoreregions"), read(file, "regionlist"),
                    read(file, "lonrange"), read(file, "latrange")
    end
end

function makeregions(regiondefinitionarray; allowmixed=false)
    regionnames, nutsdef, gadmdef = splitregiondefinitions(regiondefinitionarray)
    use_nuts, use_gadm = !all(isempty.(nutsdef)), !all(isempty.(gadmdef))
    regiontype =  (use_gadm && !use_nuts) ? :GADM  :
                  (use_nuts && !use_gadm) ? :NUTS  :
                  (use_nuts && use_gadm) ? :MIXED  :  :WEIRD
    !allowmixed && regiontype==:MIXED && error("Sorry, mixed NUTS & GADM definitions are not supported yet.")
    region = zeros(Int16, (36000,18000))    # hard code size for now
    if use_nuts
        nuts, subregionnames = read_nuts()
        makeregions_nuts!(region, nuts, subregionnames, nutsdef)
    end
    if use_gadm
        gadm, subregionnames = read_gadm()
        makeregions_gadm!(region, gadm, subregionnames, gadmdef)
    end
    return region, regiontype
end

function regions2matlab(gisregion)
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    matopen(in_datafolder("regions_$gisregion.mat"), "w", compress=true) do file
        write(file, "regions", regions)
        write(file, "offshoreregions", offshoreregions)
        write(file, "regionlist", string.(regionlist))
        write(file, "lonrange", collect(lonrange))
        write(file, "latrange", collect(latrange))
    end
end

function splitregiondefinitions(regiondefinitionarray)
    regionnames = regiondefinitionarray[:,1]
    regiondefinitions = [regdef isa Tuple ? regdef : (regdef,) for regdef in regiondefinitionarray[:,2]]
    nutsdef = [Tuple(rd for rd in regdef if rd isa NUTS) for regdef in regiondefinitions]
    gadmdef = [Tuple(rd for rd in regdef if rd isa GADM) for regdef in regiondefinitions]
    return regionnames, nutsdef, gadmdef
end

function makeregions_gadm!(region, gadm, subregionnames, regiondefinitions)
    println("Making region index matrix...")
    regionlookup = build_inverseregionlookup(regiondefinitions)
    rows, cols = size(region)
    updateprogress = Progress(cols, 1)
    for c in randperm(cols)
        for r = 1:rows
            gadm_uid = gadm[r,c]
            (gadm_uid == 0 || gadm_uid == 78413 || region[r,c] > 0) && continue    # ignore Caspian Sea (weirdly classified as a region in GADM)
            reg0, reg1, reg2 = subregionnames[gadm_uid,:]
            regid = lookup_regionnames(regionlookup, reg0, reg1, reg2)
            if regid > 0
                region[r,c] = regid
            end
        end
        next!(updateprogress)
    end
end

function makeregions_nuts!(region, nuts, subregionnames, regiondefinitions)
    println("Making region index matrix...")
    regionlookup = Dict(r => i for (i,tuptup) in enumerate(regiondefinitions)
                                    for ntup in tuptup for r in ntup.subregionnames)
    rows, cols = size(region)
    updateprogress = Progress(cols, 1)
    for c in randperm(cols)
        for r = 1:rows
            nuts_id = nuts[r,c]
            (nuts_id == 0 || region[r,c] > 0) && continue
            reg = subregionnames[nuts_id]
            while length(reg) >= 2
                regid = get(regionlookup, reg, 0)
                if regid > 0
                    region[r,c] = regid
                    break
                end
                reg = reg[1:end-1]
            end
        end
        next!(updateprogress)
    end
end

function lookup_regionnames(regionlookup, reg0, reg1, reg2)
    v = get(regionlookup, (reg0, "*", "*"), 0)
    v > 0 && return v
    v = get(regionlookup, (reg0, reg1, "*"), 0)
    v > 0 && return v
    return get(regionlookup, (reg0, reg1, reg2), 0)
end

function build_inverseregionlookup(regiondefinitions)
    d = Dict{Tuple{String,String,String}, Int}()
    for reg = 1:length(regiondefinitions)
        for regdef in regiondefinitions[reg]
            parentregions, subregionnames = regdef.parentregions, regdef.subregionnames
            regions = ["*", "*", "*"]
            regions[1:length(parentregions)] = parentregions
            for s in subregionnames
                regions[length(parentregions)+1] = s
                d[regions...] = reg
            end
        end
    end
    return d
end

function getsubregions(regtype::Type{GADM}, regionnames)
    gadm = readdlm(in_datafolder("gadmfields.csv"), ',', skipstart=1)[:, 2:4]
    for reg in regionnames
        gadm = gadm[gadm[:,1].==reg, 2:end]
    end
    return sort(unique(gadm[:,1]))
end

function getsubregions(regtype::Type{NUTS}, regionname)
    length(regionname) > 1 && error("Give only one string argument to match beginning of NUTS region names.")
    nuts = readdlm(in_datafolder("nutsfields.csv"), ',', skipstart=1)[:, 3]
    if isempty(regionname)
        return sort(unique(first.(nuts, 2)))
    else
        return sort(nuts[startswith.(nuts, regionname)])
    end
end

function subregions(regtype::Type{T} where T <: RegionType, regionnames::String...)
    reglist = join(regionnames, ", ")
    selected = string.(getsubregions(regtype, regionnames))
    selectedlist = join(selected, ", ")
    # isempty(regionnames) && println("Showing top level $regtype regions:")
    # println("$regtype($reglist): $selectedlist")
    return selected
end

function all_gadm_subregions(country::AbstractString, level::Int)
    level == 1 && return [country GADM(country)]
    return all_gadm_subregions(country, subregions(GADM, country), level)
end

function all_gadm_subregions(country::AbstractString, reg2::AbstractString, level::Int)
    reg2 == "" && return [country GADM(country)]
    level == 2 && return [reg2 GADM([country], reg2)]
    return vcat(all_gadm_subregions.(country, reg2, subregions(GADM, country, reg2))...)
end

function all_gadm_subregions(country::AbstractString, reg2::AbstractString, reg3::AbstractString)
    reg3 == "" && return [reg2 GADM([country], reg2)]
    return [reg3 GADM([country, reg2], reg3)]
end

all_gadm_subregions(countries::AbstractArray, level::Int) =
        vcat(all_gadm_subregions.(countries, level)...)

all_gadm_subregions(country::AbstractString, regions2::AbstractArray, level::Int) =
        vcat(all_gadm_subregions.(country, regions2, level)...)

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 6inch)

iris = dataset("datasets", "iris")
sp = unique(iris[:Species])
Dhl = DataFrame(yint=[3.0, 4.0, 2.5, 3.5, 2.5, 4.0], Species=repeat(sp, inner=[2]) )

function plotxg(scalexg::Scale.DiscreteScale)
    plot(iris, xgroup=:Species,
        Geom.subplot_grid(
            layer( x=:SepalLength, y=:SepalWidth, Geom.point),
            layer(Dhl, xgroup=:Species, yintercept=:yint, Geom.hline(color="red", style=:dot) ),
            ),
    scalexg,
    Theme(plot_padding=[0mm])
)
end

scales =    [Scale.xgroup(), Scale.xgroup(levels=["virginica","setosa","versicolor"]), Scale.xgroup(order=[3,1,2])]
plots = plotxg.(scales)

vstack(plots...)


text:
export solveMAP

"""
function solveMAP

Find MAP estimate using Block Coordinate Descent

Input:

	1. param       - StrainReconParam

Keyword arguments:

	2. maxTry      - number of BCD trials for solving min |Mw-d|, each trial start with random w drawn from uniform distribution
	3. storeInterm - return (M,w) with lowest misfit and history of each BCD iterations for (M,w) (storeInterm = 0)
                     return all (M,w) and intermediate history of (M,w) during BCD iteration (1)
                     return all (M,w) but no intermediate history of (M,w) during BCD iteration (2)
	4. wtolBCD     - tolerance for breaking out of BCD iteration when |w_prev - w_current| < wtolBCD
	5. maxIterBCD  - maximum number of iterations per BCD trial
	6. maxIterQP   - maximum number of iterations for findw (used in asqp)
	7. pwr         - exponent of the likelihood l_p norm

Output:

    1. Mmap - optimal M (storeInterm = 0)
              all M and intermediate history from all BCD trials, sorted by misfit (1)
              all M but no intermediate history from all BCD trials, sorted by misfit (2)
	2. wmap - optimal w (storeInterm = 0)
              all w and intermediate history from all BCD trials, sorted by misfit (1)
              all w but no intermediate history from all BCD trials, sorted by misfit (2)
	3. his  - history of misfit, |M_old-M_current|, |w_old-w_current| for optimal(M,w) with lowest misfit (0)
              history of misfit, |M_old-M_current|, |w_old-w_current| for all (M,w) and intermediate history of each BCD iterations from all BCD trials sorted by misfit (1)
              history of misfit, |M_old-M_current|, |w_old-w_current| for all (M,w) and no intermediate history from all BCD trials sorted by misfit (2)

"""

function solveMAP(param::StrainReconParam;maxTry::Int64=10,storeInterm = 0,wtolBCD::Float64=1e-3,maxIterBCD::Int64=10,maxIterQP::Int64=5,pwr::Float64=param.pwr)

	n = param.n

    if storeInterm != 1 # return (M,w) from BCD iterations with lowest misfit
        wsqp = zeros(n,maxTry)
        Msqp = falses(length(param.d),n,maxTry)
        Jsqp = zeros(maxTry)
        his = zeros(maxIterBCD,3,maxTry) # each frame for one iteration of BCD and is maxIterBCD-by-3, total of maxTry frames

        # parallelized version
        # println("parallelized version")
        updateRes = (w,m,hisf,idx) -> (wsqp[:,idx] = w;  Msqp[:,:,idx] = m; Jsqp[idx] = hisf[size(hisf,1),1]; his[1:size(hisf,1),1:size(hisf,2),idx] = hisf)

    else
        # return history of binary matrix and weights from all BCD iterations
        wsqp = zeros(n,maxIterBCD + 1,maxTry) # n-by-intermediate history for each BCD iteration-by-total BCD trials
        Msqp = falses(length(param.d),n, maxIterBCD, maxTry) # m-by-n-by-intermediate history for each BCD iteration-by-total BCD trials
        Jsqp = zeros(maxTry)
        his = zeros(maxIterBCD,3,maxTry) # each frame for one iteration of bcd-by-3, total of maxTry frames

        updateRes = (w,m,hisf,idx) -> (wsqp[1:size(w,1),1:size(w,2),idx] = w;  Msqp[1:size(m,1),1:size(m,2),1:size(m,3),idx] = m; Jsqp[idx] = hisf[size(hisf,1),1]; his[1:size(hisf,1),1:size(hisf,2),idx] = hisf)
    end

    w0   = randomw(n,maxTry)
    i=1
    nextidx = () -> (idx=i; i+=1; idx)

    @sync begin
        for pworkers = workers()
            @async begin
                while true
                    idx = nextidx()
 	 				if isempty(idx) || idx > maxTry
 	 					break
 	 				end
                    # n = 1 case, a0 is always [1]
                    if n == 1
                        m,w,hisf= remotecall_fetch(bcd,pworkers,param,[w0[idx]];wtol=wtolBCD,maxIter=maxIterBCD,maxIterQP=maxIterQP,storeInterm=storeInterm,out=-1,pwr=pwr)
                    else
                        m,w,hisf= remotecall_fetch(bcd,pworkers,param,w0[:,idx];wtol=wtolBCD,maxIter=maxIterBCD,maxIterQP=maxIterQP,storeInterm=storeInterm,out=-1,pwr=pwr)
                    end
                    updateRes(w,m,hisf,idx)
                end
            end
        end
    end

"""
    for k=1:maxTry
	    # asqp,his  = sqp(param, a0[:,k],storeInterm=false,out=-1,atol=atolSQP,maxIter=maxIterSQP, maxIterIPQP=maxIterIPQP)
	    asqp,msqp,his  = bcd(param, a0[:,k],storeInterm=false,out=-1,atol=atolSQP,maxIter=maxIterSQP,pwr=pwr)
	    Asqp[:,k] = asqp
		Msqp[:,:,k]  = msqp
	    Jsqp[k]   = misfit(msqp,asqp,param,pwr=pwr)[1];
	end

"""

	if storeInterm == 1
		id   = sortperm(Jsqp) # vector to sort by smallest to largest
		wmap = wsqp[:,:,id] # n-by-history of each BCD iteration-by-total BCD iterations
		Mmap = Msqp[:,:,:,id] # m-by-n-by-history of each BCD iteration-by-total BCD iterations, needs to be converted to bit array for input to misfit
		his = his[:,:,id] # maxIter-by-3-by-maxTry
		return Mmap,wmap,his # return all (M,w) and intermediate (M,w) during BCD iterations and history for each BCD trial
	elseif storeInterm == 0
		ksqp = indmin(Jsqp)
		wmap = wsqp[:,ksqp]
		Mmap = Msqp[:,:,ksqp] # needs to be converted to bit array for input to misfit
		his = his[:,:,ksqp]
		return Mmap,wmap,his # return (M,w) with lowest misfit and history for each BCD trial
	else
		ksqp = sortperm(Jsqp)
		wmap = wsqp[:,ksqp]
		Mmap = Msqp[:,:,ksqp] # needs to be converted to bit array for input to misfit
		his = his[:,:,ksqp]
		return Mmap,wmap,his # return all (M,w) and history for each BCD trial
	end

end