text:
"""
Visualizer for the SLAM.

Can be used to play in "live" mode as the SLAM is working
or to replay from ReplaySaver.

# Usage:

To use it in "live" preview mode, pass it as a keyword argument,
when creating SlamManager:

```julia
visualizer = Visualizer(;resolution=(900, 600), image_resolution=(370, 1226))
SlamManager(params, camera; visualizer)
```
"""
mutable struct Visualizer
    figure::Figure

    top_grid::GridLayout
    bottom_grid::GridLayout

    map_axis::Axis3
    image_axis::Makie.Axis

    positions::Observable{Vector{Point3f0}}
    position_ids::Dict{Int64, Int64}
    positions_queue::OrderedDict{Int64, Point3f0}

    image::Observable{Matrix{Gray{Float64}}}

    positions_lock::ReentrantLock
end

"""
```julia
Visualizer(;resolution, image_resolution)
```

# Arguments:

- `resolution`: Initial resolution of the figure in `(width, height)` format.
- `image_resolution`: Initial resolution of the figure in `(height, width)` format.
"""
function Visualizer(;resolution, image_resolution)
    figure = Figure(;resolution)
    top_grid = figure[1, 1:2] = GridLayout()
    bottom_grid = figure[2, 1:2] = GridLayout()

    map_axis = Axis3(
        top_grid[1, 1]; aspect=:data, azimuth=-π / 2, elevation=1)
    image_axis = Makie.Axis(
        top_grid[1, 2]; aspect=DataAspect(),
        leftspinevisible=false, rightspinevisible=false,
        bottomspinevisible=false, topspinevisible=false)

    top_grid[1, 1, Top()] = Label(figure, "Map")
    top_grid[1, 2, Top()] = Label(figure, "Current Image")

    camera_dir_element = MarkerElement(;color=:red, marker="-")
    mappoint_element = MarkerElement(;color=:gray, marker="•")

    Legend(
        bottom_grid[1, 1],
        [mappoint_element, camera_dir_element],
        ["Mappoint", "Camera position"];
        orientation=:horizontal, tellheight=true)
    trim!(figure.layout)

    positions = Observable(Point3f0[])
    position_ids = Dict{Int64, Int64}()
    positions_queue = OrderedDict{Int64, Point3f0}()

    lines!(map_axis, positions; color=:red, quality=1, linewidth=2)

    image = Observable(zeros(Gray{Float64}, image_resolution))
    image!(image_axis, image)

    Visualizer(
        figure, top_grid, bottom_grid,
        map_axis, image_axis,
        positions, position_ids, positions_queue,
        image,
        ReentrantLock())
end

Base.display(v::Visualizer) = Base.display(v.figure)

"""
```julia
set_image!(v::Visualizer, image)
```

Update image in the visualizer.
It should have the same dimensions as the original.
"""
function set_image!(v::Visualizer, image)
    v.image[] = copy!(v.image[], image)
end

"""
```julia
set_position!(v::Visualizer, position)
```

Add new position to the camera positions. This updates the plot immediately.
"""
function set_position!(v::Visualizer, position)
    v.positions[] = push!(v.positions[], position)
    autolimits!(v.map_axis)
end

"""
```julia
set_frame_wc!(v::Visualizer, frame_id, wc)
```

Add new frame `wc` to the visualizer queue.
This is used when other threads are updating the visualizer.

The queue is processed in `process_frame_wc!` method.
"""
function set_frame_wc!(v::Visualizer, frame_id, wc)
    lock(v.positions_lock) do
        base_position = SVector{4, Float64}(0, 0, 0, 1)
        position = (wc * base_position)[[1, 3, 2]]
        v.positions_queue[frame_id] = position
    end
end

"""
```julia
process_frame_wc!(v::Visualizer)
```

Process pose queue. This updates the plot immediately.
"""
function process_frame_wc!(v::Visualizer)
    lock(v.positions_lock) do
        is_done = isempty(v.positions_queue)
        do_update = !is_done

        while !is_done
            frame_id, position = popfirst!(v.positions_queue)
            is_done = isempty(v.positions_queue)

            pid = get(v.position_ids, frame_id, -1)
            if pid == -1
                v.positions[] = push!(v.positions[], position)
                v.position_ids[frame_id] = length(v.positions[])
            else
                v.positions[][pid] = position
            end
        end

        if do_update
            autolimits!(v.map_axis)
            sleep(1 / 60)
        end
    end
end


text:
using Documenter, StringBuilders

makedocs(modules=[StringBuilders],
	sitename="StringBuilders.jl",
	pages=[
        "Home" => "index.md"
    ])

deploydocs(repo="github.com/davidanthoff/StringBuilders.jl.git")

ext: jl


text:
mutable struct Box{T}
    x::T

    Box{T}() where {T} = new()
    Box{T}(x) where {T} = new(x)
end
Box(x::T) where {T} = Box{T}(x)

struct Region{N,T<:AbstractFloat,R}
    x::MVector{N,T} # center
    h::MVector{N,T} # half-width
    I::Box{R}
    E::Box{R}
    axis::Box{Int8}
end
Region(::Type{R}, x::NTuple{N,T}, h::NTuple{N,T}) where {R,N,T} =
    Region{N,T,R}(MVector(x), MVector(h), Box{R}(), Box{R}(), Box(Int8(0)))
Region(x::NTuple{N,T}, h::NTuple{N,T}) where {N,T} =
    Region{N,T,Float64}(MVector(x), MVector(h), Box(0.0), Box(0.0), Box(Int8(0)))

struct Regions{R<:Region}
    v::Vector{R}
end

struct WPoints{N,T,R}
    w::Vector{R}
    p::Vector{SVector{N,T}}
end


text:
using Manopt, ManifoldsBase, Manifolds, Test
using Random
@testset "Particle Swarm" begin
    # Test the particle swarm algorithm
    A = [1.0 3.0 4.0; 3.0 -2.0 -6.0; 4.0 -6.0 5.0]
    @testset "Eucliedean Particle Swarm" begin
        M = Euclidean(3)
        F(::Euclidean, x) = (x' * A * x) / (x' * x)
        x_start = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]
        x_start2 = deepcopy(x_start)

        Random.seed!(35)
        o = particle_swarm(M, F; x0=x_start, return_options=true)
        g = get_solver_result(o)

        Random.seed!(35)
        g2 = particle_swarm(M, F; x0=x_start2, return_options=false)
        @test isequal(g, g2)

        # the cost of g and the p[i]'s are not greater after one step
        j = argmin([F(M, y) for y in x_start])
        g0 = deepcopy(x_start[j])
        @test F(M, g) <= F(M, g0)
        for i in 1:3
            @test F(M, o.p[i]) <= F(M, x_start[i])
            # the cost of g is not greater than the cost of any p[i]
            @test F(M, g) <= F(M, o.p[i])
        end
    end
    @testset "Spherical Particle Swarm" begin
        Random.seed!(42)
        M = Sphere(2)
        F(::Sphere, x) = transpose(x) * A * x
        x_start = [random_point(M) for i in 1:3]
        v_start = [random_tangent(M, y) for y in x_start]
        p = CostProblem(M, F)
        o = ParticleSwarmOptions(x_start, v_start)
        initialize_solver!(p, o)
        step_solver!(p, o, 1)
        for i in 1:3
            # check that the new particle locations are on the manifold
            @test is_point(M, o.x[i], true)
            # check that the new velocities are tangent vectors of the original particle locations
            @test is_vector(M, o.x[i], o.velocity[i], true; atol=5 * 10^(-16))
        end
    end
end


text:
export set_raff_output_level, set_lm_output_level

"""
    residual_fg!(model::Function, data::AbstractMatrix{T}, θ::AbstractVector{T},
                 ind, r::AbstractVector{T}, rJ::AbstractMatrix{T}) where T

Compute F_I (the residual function) in `r`, using `model` and `data`,
and J_I (the Jacobian of the residual function), using `gmodel!` and
`data` in `rJ`, where I is given by `ind`, the indices associated with
the current f_I.

"""
function residual_fg!(model::Function, gmodel!::Function, data::AbstractMatrix{T}, θ::AbstractVector{T},
                      ind, r::AbstractVector{T}, rJ::AbstractMatrix{T}) where T

    for (k, i) in enumerate(ind)

        x = @view(data[i, 1:(end - 1)])

        r[k] = model(x, θ) - data[i, end]

        v = @view(rJ[k, :])

        gmodel!(v, x, θ)

    end

end

"""

    voting_strategy(model::Function, data::Array{Float64, 2}, sols::Vector{RAFFOutput}, pliminf::Int,
                    plimsup::Int)

Utility function to compute the matrix representing the voting system
used by RAFF.

It first applies a filtering strategy, to eliminate obvious local
minima, then it calculates a *magic threshold* and constructs the
distance matrix. The vector `sols` contains the solutions `s_p`, for
`p = pliminf, ... plimsup`.

"""
function voting_strategy(model::Function, data::Array{Float64, 2}, sols::Vector{RAFFOutput}, pliminf::Int,
                         plimsup::Int)

    # Remove possible stationary points, i.e., points with lower
    # values for 'p' and higher 'f'.

    eliminate_local_min!(model, data, sols)

    # Voting strategy

    lv = plimsup - pliminf + 1

    dvector = zeros(Int(lv * (lv - 1) / 2))
    dmatrix = zeros(lv, lv)
    pos = 0
    n_conv = 0

    for j = 1:lv

        # Count how many have successfully converged
        (sols[j].status == 1) && (n_conv += 1)

        for i = j + 1:lv

            dmatrix[i, j] = Inf

            if sols[i].status == 1 && sols[j].status == 1

                dmatrix[i, j] = norm(sols[i].solution - sols[j].solution)

                pos += 1

                dvector[pos] = dmatrix[i, j]

            end

        end

    end

    threshold = Inf

    if pos > 0

        dvv = @view dvector[1:pos]

        threshold = minimum(dvv) + mean(dvv) / (1.0 + sqrt(plimsup))

    elseif n_conv == 0

        @warn("No convergence for any 'p'. Returning largest.")

    end

    votsis = zeros(lv)

    @debug("Threshold: $(threshold)");

    # Actual votation

    for j = 1:lv
        # Count +1 if converged
        (sols[j].status == 1) && (votsis[j] += 1)
        # Check other distances
        for i = j + 1:lv
            if dmatrix[i, j] <=  threshold
                votsis[j] += 1
                votsis[i] += 1
            end
        end
    end

    @debug("Voting vector:", votsis)
    @debug("Distance matrix:", dmatrix)

    return votsis

end

"""

    check_ftrusted(ftrusted::Union{Float64, Tuple{Float64, Float64}}, np::Int)

Utility function to check `ftrusted` parameter in [`raff`](@ref) and
[`praff`](@ref). Throws an `ErrorException` if the percentage of
trusted points is incorrect.

"""
function check_ftrusted(ftrusted::Union{Float64, Tuple{Float64, Float64}}, np::Int)

    if typeof(ftrusted) == Float64

        (!(0.0 <= ftrusted <= 1.0)) && error("Bad value for `ftrusted`: $(ftrusted).")

        return Int(round(ftrusted * np)), np

    end

    (!(0.0 <= ftrusted[1] <= ftrusted[2] <= 1.0)) &&
        error("Bad interval for `ftrusted`: $(ftrusted[1]) > $(ftrusted[2]).")

    return Int(round(ftrusted[1] * np)), Int(round(ftrusted[2] * np))

end

"""

    eliminate_local_min!(sols::Vector{RAFFOutput})

Check if the function value of the solution found by smaller values of
`p` is not greater when compared with larger ones. This certainly
indicates that a local minimizer was found by the smaller `p`.

"""
function eliminate_local_min!(model::Function, data::Array{Float64, 2},
                              sols::Vector{RAFFOutput})

    lv = length(sols)

    sec_sol_ind = -1

    # Start from the largest p
    for i = lv:-1:1

        (sols[i].status != 1) && continue

        for j = i - 1:-1:1

            if (sols[j].status == 1) && (sols[j].f > sols[i].f)

                @debug("  Possible local minimizer for p = $(sols[j].p) " *
                       "with f = $(sols[j].f). Removing it.")

                sols[j] = RAFFOutput(0)

            end

        end

    end

    # Test the maximum 'p'

    nump, = size(data)

    maxp = sols[lv]

    if (lv > 1) && (maxp.status == 1)

        i = lv - 1

        # while (i > 0) && (sols[i].status != 1)

        #     i -= 1

        # end

        bestf = maxp.f
        j = 0

        while (i > 0)

            if  (sols[i].status == 1) && (sols[i].f < bestf)

                bestf = sols[i].f

                j = i

            end

            i -= 1

        end

        if j == 0

            @debug("  Nobody to reject p = $(maxp.p).")

        else

            sec_sol = sols[j]

            @debug("  p = $(sec_sol.p) will try to reject p = $(maxp.p).")

            nmin = 0

            for i = 1:nump

                y  = data[i, end]
                x  = @view data[i, 1:(end - 1)]
                y1 = model(x, maxp.solution)
                y2 = model(x, sec_sol.solution)

                (abs(y - y1) < abs(y - y2)) && (nmin += 1)

            end

            if nmin < lv / 2

                @debug("  nmin = $(nmin). Rejecting p = $(maxp.p).")

                sols[lv] = RAFFOutput(0)

            end

        end

    end

end

"""

This function is an auxiliary function. It finds the `p` smallest
values of vector `V` and brings them to the first `p` positions. The
indexes associated with the `p` smallest values are stored in `ind`.

"""
function sort_fun!(V::Vector{Float64}, ind::Vector{Int}, p::Int)

    # If p is invalid, returns an empty view
    (p <= 0) && (return @view(ind[1:p]), @view(V[1:p]))

    npun = length(ind)

    for i = 1:npun
        ind[i] = i
    end

    for i = 1:p

        for j = i + 1:npun

            if V[i] > V[j]

                aux = ind[j]
                ind[j] = ind[i]
                ind[i] = aux

                vaux = V[j]
                V[j] = V[i]
                V[i] = vaux
            end
        end
    end

    return @view(ind[1:p]), @view(V[1:p])

end

"""

    set_raff_output_level(level::LogLevel)

Set the output level of [`raff`](@ref) and [`praff`](@ref) algorithms
to the desired logging level. Options are (from highly verbose to just
errors): `Logging.Debug`, `Logging.Info`, `Logging.Warn` and
`Logging.Error`. The package
[`Logging`](https://docs.julialang.org/en/v1.0/stdlib/Logging/index.html)
needs to be loaded.

Defaults to `Logging.Error`.

"""
function set_raff_output_level(level::LogLevel)

    global raff_logger.x = ConsoleLogger(stdout, level)

end

"""

    set_lm_output_level(level::LogLevel)

Set the output level of [`lmlovo`](@ref) algorithm to the desired
logging level. Options are (from highly verbose to just errors):
`Logging.Debug`, `Logging.Info`, `Logging.Warn` and
`Logging.Error`. The package
[`Logging`](https://docs.julialang.org/en/v1.0/stdlib/Logging/index.html)
needs to be loaded.

Defaults to `Logging.Error`.

"""
function set_lm_output_level(level::LogLevel)

    global lm_logger.x = ConsoleLogger(stdout, level)

end


text:
using .DimensionalData: DimensionalArray, DimensionalData, data, Dim, metadata

_dname(::DimensionalData.Dim{N}) where N = N
dimname(x::DimensionalArray, i) = _dname(DimensionalData.dims(x)[i])

dimvals(x::DimensionalArray,i) = DimensionalData.dims(x)[i].val

getdata(x::DimensionalArray) = data(x)

getattributes(x::DimensionalArray) = metadata(x)

function yaxcreate(::Type{<:DimensionalArray},data,dnames,dvals,atts)
  d = ntuple(ndims(data)) do i
    Dim{Symbol(dnames[i])}(dvals[i])
  end
  DimensionalArray(data,d,metadata = atts)
end


text:
# Autogenerated wrapper script for CImPlot_jll for x86_64-w64-mingw32
export compile_commands, libcimplot

using CImGui_jll
JLLWrappers.@generate_wrapper_header("CImPlot")
JLLWrappers.@declare_file_product(compile_commands)
JLLWrappers.@declare_library_product(libcimplot, "libcimplot.dll")
function __init__()
    JLLWrappers.@generate_init_header(CImGui_jll)
    JLLWrappers.@init_file_product(
        compile_commands,
        "share\\compile_commands.json",
    )

    JLLWrappers.@init_library_product(
        libcimplot,
        "bin\\libcimplot.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const EXTRA_KW_MAKE = """

* `eltype=Float64`: machine type of points (any subtype of
   `AbstractFloat`).

* `rng=nothing`: any `AbstractRNG` object, or integer to seed a
  `MersenneTwister` (for reproducibility).

* `as_table=true`: whether to return the points as a table (true)
  or a matrix (false). """

const EXTRA_CLASSIFICATION =
    "If `false` the target `y` has integer element type. "

"""
    finalize_Xy(X, y, shuffle, as_table, eltype, rng; clf)

Internal function to  finalize the `make_*` functions.

"""
function finalize_Xy(X, y, shuffle, as_table, eltype, rng; clf::Bool=true)
    # Shuffle the rows if required
    if shuffle
                X, y = shuffle_rows(X, y; rng=rng)
        end
        if eltype != Float64
                X = convert.(eltype, X)
        end
        # return as matrix if as_table=false
        as_table || return X, y
        clf && return MLJBase.table(X), categorical(y)
        return MLJBase.table(X), y
end

### CLASSIFICATION TOY DATASETS

"""
    runif_ab(rng, n, p, a, b)

Internal function to generate `n` points in `[a, b]ᵖ` uniformly at random.

"""
runif_ab(rng, n, p, a, b) = (b - a) .* rand(rng, n, p) .+ a

"""
    X, y = make_blobs(n=100, p=2; kwargs...)

Generate Gaussian blobs for clustering and classification
problems.

### Return value

By default, a table `X` with `p` columns (features) and `n` rows
(observations), together with a corresponding vector of `n`
`Multiclass` target observations `y`, indicating blob membership.

### Keyword arguments

* `shuffle=true`: whether to shuffle the resulting points,

* `centers=3`: either a number of centers or a `c x p` matrix with `c`
  pre-determined centers,

* `cluster_std=1.0`: the standard deviation(s) of each blob,

* `center_box=(-10. => 10.)`: the limits of the `p`-dimensional cube
  within which the cluster centers are drawn if they are not provided,
  $(EXTRA_KW_MAKE*EXTRA_CLASSIFICATION)

### Example

```
X, y = make_blobs(100, 3; centers=2, cluster_std=[1.0, 3.0])
```

"""
function make_blobs(n::Integer=100,
                    p::Integer=2;
                    shuffle::Bool=true,
                    centers::Union{<:Integer,Matrix{<:Real}}=3,
                    cluster_std::Union{<:Real,Vector{<:Real}}=1.0,
                    center_box::Pair{<:Real,<:Real}=(-10.0 => 10.0),
                    as_table::Bool=true,
                    eltype::Type{<:AbstractFloat}=Float64,
                    rng=nothing)

    # check arguments make sense
    if n < 1 || p < 1
        throw(ArgumentError(
            "Expected `n` and `p` to be at least 1."))
    end
    if center_box.first >= center_box.second
        throw(ArgumentError(
            "Domain for the centers improperly defined expected a pair " *
            "`a => b` with `a < b`."))
    end

    if rng === nothing
        rng = Random.GLOBAL_RNG
    elseif rng isa Integer
        rng = Random.MersenneTwister(rng)
    end

    if centers isa Matrix
        if size(centers, 2) != p
            throw(ArgumentError(
                "The centers provided have dimension ($(size(centers, 2))) " *
                "that doesn't match the one specified ($(p))."))
        end
        n_centers = size(centers, 1)
    else
        # in case the centers aren't provided, draw them from the box
        n_centers = centers
        centers   = runif_ab(rng, n_centers, p, center_box...)
    end
    if cluster_std isa Vector
        if length(cluster_std) != n_centers
            throw(ArgumentError(
                "$(length(cluster_std)) standard deviations given but there " *
                "are $(n_centers) centers."))
        end
        if any(cluster_std .<= 0)
            throw(ArgumentError(
                "Cluster(s) standard deviation(s) must be positive."))
        end
    else
        # In case only one std is given, repeat it for each center
        cluster_std = fill(cluster_std, n_centers)
    end

    # split points equally among centers
    ni, r    = divrem(n, n_centers)
    ns       = fill(ni, n_centers)
    ns[end] += r

    # vector of memberships
    y = vcat((fill(i, ni) for (i, ni) in enumerate(ns))...)

    # Pre-generate random points then modify for each center
    X    = randn(rng, n, p)
    nss  = [1, cumsum(ns)...]
    # ranges of rows for each center
    rows = [nss[i]:nss[i+1] for i in 1:n_centers]
    @inbounds for c in 1:n_centers
        Xc = view(X, rows[c], :)
        # adjust standard deviation
        Xc .*= cluster_std[c]
        # adjust center
        Xc .+= centers[c, :]'
    end

    return finalize_Xy(X, y, shuffle, as_table, eltype, rng)
end

"""
    X, y = make_circles(n=100; kwargs...)

Generate `n` labeled points close to two concentric circles for
classification and clustering models.

### Return value

By default, a table `X` with `2` columns and `n` rows (observations),
together with a corresponding vector of `n` `Multiclass` target
observations `y`. The target is either `0` or `1`, corresponding to
membership to the smaller or larger circle, respectively.

### Keyword arguments

* `shuffle=true`: whether to shuffle the resulting points,

* `noise=0`: standard deviation of the Gaussian noise added to the data,

* `factor=0.8`: ratio of the smaller radius over the larger one,

$(EXTRA_KW_MAKE*EXTRA_CLASSIFICATION)

### Example

```
X, y = make_circles(100; noise=0.5, factor=0.3)
```

"""
function make_circles(n::Integer=100;
                      shuffle::Bool=true,
                      noise::Real=0.,
                      factor::Real=0.8,
                      as_table::Bool=true,
                      eltype::Type{<:AbstractFloat}=Float64,
                      rng=nothing)

    # check arguments make sense
    if n < 1
        throw(ArgumentError(
            "Expected `n` to be at least 1."))
    end
    if noise < 0
        throw(ArgumentError(
            "Noise argument cannot be negative."))
    end
    if !(0 < factor < 1)
        throw(ArgumentError(
            "Factor argument must be strictly between 0 and 1."))
    end

        if rng === nothing
                rng = Random.GLOBAL_RNG
        elseif rng isa Integer
                rng = Random.MersenneTwister(rng)
        end

    # Generate points on a 2D circle
    θs = runif_ab(rng, n, 1, 0, 2pi)

    n0 = div(n, 2)

    X = hcat(cos.(θs), sin.(θs))
    X[1:n0, :] .*= factor

    y = ones(Int, n)
    y[1:n0] .= 0

    if !iszero(noise)
        X .+= noise .* randn(rng, n, 2)
    end

        return finalize_Xy(X, y, shuffle, as_table, eltype, rng)
end

"""
        make_moons(n::Int=100; kwargs...)

Generates labeled two-dimensional points lying close to two
interleaved semi-circles, for use with classification and clustering
models.

### Return value

By default, a table `X` with `2` columns and `n` rows (observations),
together with a corresponding vector of `n` `Multiclass` target
observations `y`. The target is either `0` or `1`, corresponding to
membership to the left or right semi-circle.

### Keyword arguments

* `shuffle=true`: whether to shuffle the resulting points,

* `noise=0.1`: standard deviation of the Gaussian noise added to the data,

* `xshift=1.0`: horizontal translation of the second center with respect to
  the first one.

* `yshift=0.3`: vertical translation of the second center with respect
  to the first one.  $(EXTRA_KW_MAKE*EXTRA_CLASSIFICATION)

### Example

```
X, y = make_moons(100; noise=0.5)
```

"""
function make_moons(n::Int=150;
                    shuffle::Bool=true,
                    noise::Real=0.1,
                    xshift::Real=1.0,
                    yshift::Real=0.3,
                    as_table::Bool=true,
                    eltype::Type{<:AbstractFloat}=Float64,
                    rng=nothing)

    # check arguments make sense
    if n < 1
        throw(ArgumentError(
            "Expected `n` to be at least 1."))
    end
    if noise < 0
        throw(ArgumentError(
            "Noise argument cannot be negative."))
    end

        if rng === nothing
                rng = Random.GLOBAL_RNG
        elseif rng isa Integer
                rng = Random.MersenneTwister(rng)
        end

    n1 = div(n, 2)
    n2 = n - n1

    θs = runif_ab(rng, n, 1, 0, pi)
    θs[n2+1:end] .*= -1

    X = hcat(cos.(θs), sin.(θs))

    X[n2+1:end, 1] .+= xshift
    X[n2+1:end, 2] .+= yshift

    y = ones(Int, n)
    y[1:n1] .= 0

    if !iszero(noise)
        X .+= noise .* randn(rng, n, 2)
    end

    return finalize_Xy(X, y, shuffle, as_table, eltype, rng)
end

### REGRESSION TOY DATASETS

"""
    augment_X(X, fit_intercept)

Given a matrix `X`, append a column of ones if `fit_intercept` is true.
See [`make_regression`](@ref).

"""
function augment_X(X::Matrix{<:Real}, fit_intercept::Bool)
        fit_intercept || return X
        return hcat(X, ones(eltype(X), size(X, 1)))
end

"""
    sparsify!(rng, θ, s)

Make portion `s` of vector `θ` exactly 0.

"""
sparsify!(rng, θ, s) =
        (θ .*= (rand(rng, length(θ)) .< s))

"""Add outliers to portion s of vector."""
outlify!(rng, y, s) =
        (n = length(y); y .+= 20 * randn(rng, n) .* (rand(rng, n) .< s))

const SIGMOID_64 = log(Float64(1)/eps(Float64) - Float64(1))
const SIGMOID_32 = log(Float32(1)/eps(Float32) - Float32(1))

"""
    sigmoid(x)

Return the sigmoid computed in a numerically stable way:

``σ(x) = 1/(1+exp(-x))``

"""
function sigmoid(x::Float64)
        x > SIGMOID_64  && return one(x)
        x < -SIGMOID_64 && return zero(x)
        return one(x) / (one(x) + exp(-x))
end
function sigmoid(x::Float32)
        x > SIGMOID_32  && return one(x)
        x < -SIGMOID_32 && return zero(x)
        return one(x) / (one(x) + exp(-x))
end
sigmoid(x) = sigmoid(float(x))

"""
    make_regression(n, p; kwargs...)

Generate Gaussian input features and a linear response with Gaussian
noise, for use with regression models.

### Return value

By default, a table `X` with `p` columns and `n` rows (observations),
together with a corresponding vector of `n` `Continuous` target
observations `y`.

### Keywords

* `intercept=true: whether to generate data from a model with
  intercept,

* `sparse=0`: portion of the generating weight vector that is sparse,

* `noise=0.1`: standard deviation of the Gaussian noise added to the
  response,

* `outliers=0`: portion of the response vector to make as outliers by
  adding a random quantity with high variance. (Only applied if
  `binary` is `false`)

* `binary=false`: whether the target should be binarized (via a sigmoid).
$EXTRA_KW_MAKE

### Example

```
X, y = make_regression(100, 5; noise=0.5, sparse=0.2, outliers=0.1)
```

"""
function make_regression(n::Int=100,
                         p::Int=2;
                         intercept::Bool=true,
                         sparse::Real=0,
                         noise::Real=0.1,
                         outliers::Real=0,
                         binary::Bool=false,
                         as_table::Bool=true,
                         eltype::Type{<:AbstractFloat}=Float64,
                         rng=nothing)

    # check arguments make sense
    if n < 1 || p < 1
        throw(ArgumentError(
            "Expected `n` and `p` to be at least 1."))
    end
    if !(0 <= sparse < 1)
        throw(ArgumentError(
            "Sparsity argument must be in [0, 1)."))
    end
    if noise < 0
        throw(ArgumentError(
            "Noise argument cannot be negative."))
    end
    if !(0 <= outliers <= 1)
        throw(ArgumentError(
            "Outliers argument must be in [0, 1]."))
    end

    if rng === nothing
        rng = Random.GLOBAL_RNG
    elseif rng isa Integer
        rng = Random.MersenneTwister(rng)
    end

    X = augment_X(randn(rng, n, p), intercept)
    θ = randn(rng, p + Int(intercept))
    sparse > 0 && sparsify!(rng, θ, sparse)
    y = X * θ

    if !iszero(noise)
        y .+= noise .* randn(rng, n)
    end

    if binary
        y = rand(rng, n) .< sigmoid.(y)
    else
        if !iszero(outliers)
            outlify!(rng, y, outliers)
        end
    end

    return finalize_Xy(X[:,1:end-Int(intercept)], y, false,
                       as_table, eltype, rng; clf=binary)
end


text:
module Integrate

import DifferentialEquations: solve, ODEProblem, Tsit5

const alg = Tsit5()

function solve(f, a::Float64, b::Float64; kwargs...)
    solve(ODEProblem((u, p, t) -> f(t), 0.0, (a, b)), alg; kwargs...)
end

"""
Integrate f from lower bound a to upper bound b.
"""
function integrate(f, a::Float64, b::Float64)::Float64
    return solve(f, a, b; save_everystep=false)[end]
end

function integrate(f, a::Float64, b::Array{Float64, 1})::Array{Float64, 1}
    bmax = maximum(b)
    return solve(f, a, bmax; saveat = b).u
end

end


text:
"""
    DecisionVariable <: AbstractVariableRef

Identifier for a decision in a stochastic program. Holds a reference to the stochastic program, the stage the decision is taken in, and its corresponding MOI.VariableIndex.
"""
struct DecisionVariable <: JuMP.AbstractVariableRef
    stochasticprogram::StochasticProgram
    stage::Int
    index::MOI.VariableIndex
end

# Getters #
# ========================== #
"""
    decision(dvar::DecisionVariable)

Return the internal `Decision` associated with the first-stage `dvar`.
"""
function decision(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `decision(dvar, scenario_index)`.")
    return decision(structure(owner_model(dvar)), index(dvar), stage(dvar))
end
"""
    decision(dvar::DecisionVariable, scenario_index::Integer)

Return the scenario-dependent internal `Decision` associated with `dvar` at `scenario_index`.
"""
function decision(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `decision(dvar)`.")
    return decision(structure(owner_model(dvar)), index(dvar), stage(dvar), scenario_index)
end
"""
    stage(dvar::DecisionVariable)

Return the stage of `dvar`.
"""
function stage(dvar::DecisionVariable)
    return dvar.stage
end
"""
    state(dvar::DecisionVariable)

Return the `DecisionState` of the first-stage `dvar`.
"""
function state(dvar::DecisionVariable)
    return decision(dvar).state
end
"""
    state(dvar::DecisionVariable, scenario_index::Integer)

Return the scenario-dependent `DecisionState` of `dvar` at `scenario_index`.
"""
function state(dvar::DecisionVariable, scenario_index::Integer)
    return decision(dvar,scenario_index).state
end

# Setters #
# ========================== #
function take_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, vals::AbstractVector)
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, stochasticprogram), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are first-stage
    stage(dvars[1]) == 1 || error("Decisions are scenario-dependent, consider `take_decisions!(sp, dvars, vals, scenario_index)`")
    # Check decision length
    length(dvars) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(dvars)).")
    # Update decisions
    for (dvar, val) in zip(dvars, vals)
        fix(dvar, val)
    end
    return nothing
end
function take_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, vals::AbstractVector, scenario_index::Integer)
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, stochasticprogram), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are not first-stage
    stage(dvars[1]) > 1 || error("Decisions are not scenario-dependent, consider `take_decisions!(sp, dvars, vals)`")
    # Check decision length
    length(dvars) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(dvars)).")
    # Update decisions
    for (dvar, val) in zip(dvars, vals)
        fix(dvar, scenario_index, val)
    end
    return nothing
end

function untake_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable})
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, model), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are first-stage
    stage(dvars[1]) == 1 || error("Decisions are scenario-dependent, consider `unttake_decisions!(sp, dvars, scenario_index)`")
    # Update decisions
    need_update = false
    for dvar in dvars
        unfix(dvar)
    end
    return nothing
end
function untake_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, scenario_index::Integer)
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, model), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are first-stage
    stage(dvars[1]) > 1 || error("Decisions are not scenario-dependent, consider `take_decisions!(sp, dvars, vals)`")
    # Update decisions
    for dvar in dvars
        unfix(dvar, scenario_index)
    end
    return nothing
end

# MOI #
# ========================== #
function MOI.get(stochasticprogram::StochasticProgram, attr::MOI.AbstractVariableAttribute,
                 dvar::DecisionVariable)
    check_belongs_to_model(dvar, stochasticprogram)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        if haskey(cache, :solution)
            # Returned cached solution if possible
            try
                return MOI.get(cache[:solution], attr, index(dvar))
            catch
            end
        end
        if haskey(cache, :node_solution_1)
            # Value was possibly only cached in first-stage solution
            try
                return MOI.get(cache[:node_solution_1], attr, index(dvar))
            catch
            end
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        return MOI.get(optimizer(stochasticprogram), attr, index(dvar))
    end
    return MOI.get(backend(proxy(stochasticprogram, stage(dvar))), attr, index(dvar))
end
function MOI.get(stochasticprogram::StochasticProgram, attr::ScenarioDependentVariableAttribute,
                 dvar::DecisionVariable)
    check_belongs_to_model(dvar, stochasticprogram)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        key = Symbol(:node_solution_, attr.stage, :_, attr.scenario_index)
        if haskey(cache, key)
            # Returned cached solution if possible
            try
                return MOI.get(cache[key], attr.attr, index(dvar))
            catch
            end
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        try
            # Try to get scenario-dependent value directly
            return MOI.get(optimizer(stochasticprogram), attr, index(dvar))
        catch
            # Fallback to resolving scenario-dependence in structure if
            # not supported natively by optimizer
            MOI.get(structure(stochasticprogram), attr, index(dvar))
        end
    end
    # Get value from structure if not set by optimizer
    return MOI.get(structure(stochasticprogram), attr, index(dvar))
end

function MOI.set(stochasticprogram::StochasticProgram, attr::MOI.AbstractVariableAttribute,
                 dvar::DecisionVariable, value)
    check_belongs_to_model(dvar, stochasticprogram)
    MOI.set(structure(stochasticprogram), attr, index(dvar), value)
    return nothing
end

# JuMP variable interface #
# ========================== #
"""
    name(dvar::DecisionVariable, scenario_index::Integer)::String

Get the name of the decision variable `dvar`.
"""
function JuMP.name(dvar::DecisionVariable)
    return MOI.get(owner_model(dvar), MOI.VariableName(), dvar)::String
end
"""
    name(dvar::DecisionVariable, scenario_index::Integer)::String

Get the name of the scenario-dependent decision variable `dvar` in scenario `scenario_index`.
"""
function JuMP.name(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `name(dvar)`.")
    attr = ScenarioDependentVariableAttribute(stage(dvar), scenario_index, MOI.VariableName())
    return MOI.get(owner_model(dvar), attr, dvar)::String::String
end
"""
    set_name(dvar::DecisionVariable, scenario_index::Integer, name::String)

Set the name of the decision variable `dvar` to `name`.
"""
function JuMP.set_name(dvar::DecisionVariable, name::String)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_name(dvar, scenario_index, name)`.")
    return MOI.set(owner_model(dvar), MOI.VariableName(), dvar, name)
end
"""
    set_name(dvar::DecisionVariable, scenario_index::Integer, name::String)

Set the name of the scenario-dependent decision variable `dvar` in scenario `scenario_index` to `name`.
"""
function JuMP.set_name(dvar::DecisionVariable, scenario_index::Integer, name::String)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_name(dvar, name)`.")
    attr = ScenarioDependentVariableAttribute(stage(dvar), scenario_index, MOI.VariableName())
    return MOI.set(owner_model(dvar), attr, dvar, name)
end
"""
    decision_by_name(stochasticprogram::Stochasticprogram,
                     stage::Integer,
                     name::String)::Union{AbstractVariableRef, Nothing}

Returns the reference of the variable with name attribute `name` at `stage` of `stochasticprogram` or `Nothing` if
no variable has this name attribute. Throws an error if several variables have
`name` as their name attribute at stage `s`.
"""
function decision_by_name(stochasticprogram::StochasticProgram{N}, stage::Integer, name::String) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    index = MOI.get(backend(proxy(stochasticprogram, stage)), MOI.VariableIndex, name)
    if index isa Nothing
        return nothing
    else
        return DecisionVariable(stochasticprogram, stage, index)
    end
end
"""
    index(dvar::DecisionVariable)::MOI.VariableIndex

Return the index of the decision variable that corresponds to `dvar` in the MOI backend.
"""
JuMP.index(dvar::DecisionVariable) = dvar.index
"""
    optimizer_index(dvar::DecisionVariable)::MOI.VariableIndex

Return the index of the variable that corresponds to `dvar` in the optimizer model.
"""
function JuMP.optimizer_index(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `optimizer_index(dvar, scenario_index)`.")
    return JuMP._moi_optimizer_index(structure(owner_model(dvar)), index(dvar))
end
"""
    optimizer_index(dvar::DecisionVariable, scenario_index)::MOI.VariableIndex

Return the index of the variable that corresponds to the scenario-dependent `dvar` in the optimizer model at `scenario_index`.
"""
function JuMP.optimizer_index(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `optimizer_index(dvar)`.")
    return JuMP._moi_optimizer_index(structure(owner_model(dvar)), index(dvar), scenario_index)
end
"""
    has_duals(stochasticprogram::StochasticProgram; result::Int = 1)

Return `true` if the solver has a primal solution in the first-stage of
`stochasticprogram` in result index `result` available to query,
otherwise return `false`.
"""
function JuMP.has_values(stochasticprogram::StochasticProgram; result::Int = 1)
    return primal_status(stochasticprogram; result = result) != MOI.NO_SOLUTION
end
"""
    has_values(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)

Return `true` if the solver has a primal solution in the node at stage
`stage` and scenario `scenario_index` in result index `result`
available to query, otherwise return `false`.
"""
function JuMP.has_values(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)
    return primal_status(stochasticprogram, stage, scenario_index; result = result) != MOI.NO_SOLUTION
end
"""
    has_values(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)

Return `true` if the solver has a primal solution in scenario `scenario_index`
in result index `result` available to query, otherwise return `false`.
"""
function has_values(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)
    return has_values(stochasticprogram, 2, scenario_index; result = result)
end
"""
    value(dvar::DecisionVariable; result = 1)

Return the value of the first-stage decision variable `dvar`
associated with result index `result` of the most-recent
returned by the solver.
"""
function JuMP.value(dvar::DecisionVariable; result::Int = 1)::Float64
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `value(dvar, scenario_index)`.")
    d = decision(dvar)
    if d.state == Taken
        # If decision has been fixed the value can be fetched
        # directly
        return d.value
    end
    return MOI.get(owner_model(dvar), MOI.VariablePrimal(result), dvar)
end
"""
    value(dvar::DecisionVariable, scenario_index::Integer; result = 1)

Return the scenario-dependent value of the decision variable `dvar`
associated with result index `result` at `scenario_index` of the
most-recent returned by the solver.
"""
function JuMP.value(dvar::DecisionVariable, scenario_index::Integer; result::Int = 1)::Float64
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_fixed(dvar)`.")
    d = decision(dvar, scenario_index)
    if d.state == Taken
        # If decision has been fixed the value can be fetched
        # directly
        return d.value
    end
    attr = ScenarioDependentVariableAttribute(stage(dvar), scenario_index, MOI.VariablePrimal(result))
    return MOI.get(owner_model(dvar), attr, dvar)
end
"""
    value(dvar_expr::Union{GenericAffExpr{T,DecisionVariable}, GenericQuadExpr{T,DecisionVariable}}, stage_to_scenario::Dict{Int,Int}) where T

Evaluate `dvar_expr` where the value of a given `dvar` is found in the scenario returned by the provided `stage_to_scenario` map.
"""
function JuMP.value(dvar_expr::Union{GenericAffExpr{T,DecisionVariable}, GenericQuadExpr{T,DecisionVariable}}, stage_to_scenario::Dict{Int,Int}; result::Int = 1)::Float64 where T
    var_value = (dvar) -> begin
        if stage(dvar) == 1
            return value(dvar; result = result)
        end
        return value(dvar, stage_to_scenario[stage(dvar)]; result = result)
    end
    return value(dvar_expr, var_value)
end
"""
    reduced_cost(dvar::DecisionVariable)::Float64

Return the reduced cost associated with the decision variable `dvar`.
"""
function JuMP.reduced_cost(dvar::DecisionVariable)::Float64
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `reduced_cost(dvar, scenario_index)`.")
    sp = owner_model(dvar)
    if !has_duals(sp)
        error("Unable to query reduced cost of variable because stochastic program does" *
              " not have duals available in the first stage.")
    end
    sign = objective_sense(sp) == MOI.MIN_SENSE ? 1.0 : -1.0
    if is_fixed(dvar)
        return sign * dual(FixRef(dvar))
    end
    rc = 0.0
    if has_upper_bound(dvar)
        rc += dual(UpperBoundRef(dvar))
    end
    if has_lower_bound(dvar)
        rc += dual(LowerBoundRef(dvar))
    end
    return sign * rc
end
"""
    reduced_cost(dvar::DecisionVariable)::Float64

Return the reduced cost associated with the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.reduced_cost(dvar::DecisionVariable, scenario_index::Integer)::Float64
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `reduced_cost(dvar)`.")
    sp = owner_model(dvar)
    if !has_duals(sp, scenario_index)
        error("Unable to query reduced cost of variable because stochastic program does" *
              " not have duals available in scenario `scenario_index`.")
    end
    sign = objective_sense(sp) == MOI.MIN_SENSE ? 1.0 : -1.0
    if is_fixed(dvar, scenario_index)
        return sign * dual(FixRef(dvar), scenario_index)
    end
    rc = 0.0
    if has_upper_bound(dvar, scenario_index)
        rc += dual(UpperBoundRef(dvar), scenario_index)
    end
    if has_lower_bound(dvar, scenario_index)
        rc += dual(LowerBoundRef(dvar), scenario_index)
    end
    return sign * rc
end
"""
    is_fixed(dvar::DecisionVariable)

Return `true` if `dvar` is a fixed first-stage decision variable.
"""
function JuMP.is_fixed(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_fixed(dvar, scenario_index)`.")
    if state(dvar) == Taken
        return true
    end
    return false
end
"""
    is_fixed(dvar::DecisionVariable, scenario_index)

Return `true` if `dvar` is a fixed decision variable in `scenario_index`.
"""
function JuMP.is_fixed(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_fixed(dvar)`.")
    if state(dvar, scenario_index) == Taken
        return true
    end
    return false
end
"""
    fix(dvar::DecisionVariable, val::Number)

Fix the first-stage decision associated with `dvar` to `val`. In contexts
where `dvar` is a variable, the variable is fixed to the value. In
contexts where `dvar` is a known parameter value, the value is updated.
"""
function JuMP.fix(dvar::DecisionVariable, val::Number)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `fix(dvar, scenario_index, val)`.")
    fix(structure(owner_model(dvar)), index(dvar), stage(dvar), val)
    return nothing
end
"""
    fix(dvar::DecisionVariable, scenario_index::Integer, val::Number)

Fix the scenario-dependent decision associated with `dvar` at `scenario_index`
to `val`. In contexts where `dvar` is a variable, the variable is fixed to the
value. In contexts where `dvar` is a known parameter value, the value is updated.
"""
function JuMP.fix(dvar::DecisionVariable, scenario_index::Integer, val::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `fix(dvar, val)`.")
    fix(structure(owner_model(dvar)), index(dvar), stage(dvar), scenario_index, val)
    return nothing
end
"""
    unfix(dvar::DecisionVariable)

Unfix the first-stage decision associated with `dvar`. If the decision
is already in a `NotTaken` state, this does nothing.
"""
function JuMP.unfix(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `unfix(dvar, scenario_index)`.")
    if state(dvar) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    unfix(structure(owner_model(dvar)), index(dvar), stage(dvar))
    return nothing
end
"""
    unfix(dvar::DecisionVariable, scenario_index::Integer)

Unfix the scenario-dependent decision associated with `dvar`
at `scenario_index`. If the decision is already in a `NotTaken`
state, this does nothing.
"""
function JuMP.unfix(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `unfix(dvar)`.")
    if state(dvar, scenario_index) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    unfix(structure(owner_model(dvar)), index(dvar), stage(dvar), scenario_index)
    return nothing
end

"""
    FixRef(dvar::DecisionVariable)

Return a constraint reference to the constraint fixing the value of the decision `dvar`. Errors
if one does not exist.
"""
function JuMP.FixRef(dvar::DecisionVariable)
    moi_fix =  MOI.ConstraintIndex{SingleDecision, MOI.EqualTo{Float64}}
    return SPConstraintRef{moi_fix, ScalarShape}(owner_model(dvar),
                                                 stage(dvar),
                                                 moi_fix(index(dvar).value),
                                                 ScalarShape())
end

JuMP.owner_model(dvar::DecisionVariable) = dvar.stochasticprogram

struct DecisionVariableNotOwned <: Exception
    dvar::DecisionVariable
end

function JuMP.check_belongs_to_model(dvar::DecisionVariable, stochasticprogram::StochasticProgram)
    if owner_model(dvar) !== stochasticprogram
        throw(DecisionNotOwned(dvar))
    end
end
"""
    is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable)

Return `true` if `dvar` refers to a valid first-stage decision variable in `stochasticprogram`.
"""
function JuMP.is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_valid(dvar, scenario_index)`.")
    return stochasticprogram === owner_model(dvar) &&
        MOI.is_valid(structure(stochasticprogram), index(dvar), stage(dvar))
end
"""
    is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent `dvar` refers to a valid decision variable in `stochasticprogram` at `scenario_index`.
"""
function JuMP.is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_valid(dvar)`.")
    return stochasticprogram === owner_model(dvar) &&
        MOI.is_valid(structure(stochasticprogram), index(dvar), stage(dvar), scenario_index)
end
"""
    delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable)

Delete the first-stage decision variable associated with `dvar` from the `stochasticprogram`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `delete(stochasticprogram, dvar, scenario_index)`.")
    if stochasticprogram !== owner_model(dvar)
        error("The decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    proxy_ = proxy(stochasticprogram, stage(dvar))
    JuMP.delete(proxy_, DecisionRef(proxy_, index(dvar)))
    MOI.delete(structure(stochasticprogram), index(dvar), stage(dvar))
    return nothing
end
"""
    delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)

Delete the scenario-dependent decision variable associated with `dvar` from the `stochasticprogram` at `scenario_index`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `delete(stochasticprogram, dvar)`.")
    if stochasticprogram !== owner_model(dvar)
        error("The decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    MOI.delete(structure(stochasticprogram), index(dvar), stage(dvar), scenario_index)
    return nothing
end
"""
    delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable})

Delete the decisions associated with `dvars` from the `stochasticprogram`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable})
    all(stage.(dvars) .== stage(dvars[1])) || error("$dvars are not all from the same stage")
    stage(dvars[1]) > 1 && error("Some of $dvars are scenario dependent, consider `delete(stochasticprogram, dvars, scenario_index)`.")
    if any(stochasticprogram !== owner_model(dvar) for dvar in dvars)
        error("A decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    proxy_ = proxy(stochasticprogram, stage(dvars[1]))
    JuMP.delete(proxy_, DecisionRef.(proxy_, index.(dvars)))
    MOI.delete(structure(stochasticprogram), index.(dvars), stage(dvars[1]))
    return nothing
end
"""
    delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, scenario_index::Integer)

Delete the scenario-dependent decisions associated with `dvars` from the `stochasticprogram` at `scenario_index`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, scenario_index::Integer)
    all(stage.(dvars) .== stage(dvars[1])) || error("$dvars are not all from the same stage")
    stage(dvars[1]) == 1 && error("$dvars are not scenario dependent, consider `delete(stochasticprogram, dvars)`.")
    if any(stochasticprogram !== owner_model(dvar) for dvar in dvars)
        error("A decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    MOI.delete(structure(stochasticprogram), index.(dvars), stage(dvars[1]), scenario_index)
    return nothing
end
"""
    has_lower_bound(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` has a lower bound.
"""
function JuMP.has_lower_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `has_lower_bound(dvar, scenario_index)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar))
end
"""
    has_lower_bound(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` has a lower bound at `scenario_index`.
"""
function JuMP.has_lower_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `has_lower_bound(dvar)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar), scenario_index)
end
"""
    LowerBoundRef(dvar::DecisionVariable)

Return a constraint reference to the lower bound constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.LowerBoundRef(dvar::DecisionVariable)
    moi_lb =  MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}
    return SPConstraintRef{moi_lb, ScalarShape}(owner_model(dvar),
                                                stage(dvar),
                                                moi_lb(index(dvar).value),
                                                ScalarShape())
end
"""

    set_lower_bound(dvar::DecisionVariable)

Set the lower bound of the first-stage decision variable `dvar` to `lower`. If one does not exist, create a new lower bound constraint.
"""
function JuMP.set_lower_bound(dvar::DecisionVariable, lower::Number)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_lower_bound(dvar, scenario_index, lower)`.")
    new_set = MOI.GreaterThan(convert(Float64, lower))
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if has_lower_bound(dvar)
        # Update existing bound
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dvar).value)
        MOI.set(structure(owner_model(dvar)), MOI.ConstraintSet(), cindex, new_set)
        # Update proxy
        MOI.set(backend(proxy_), MOI.ConstraintSet(), cindex, new_set)
    else
        # Add new lower bound constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), new_set)
        # Update proxy
        MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), new_set)
    end
    return nothing
end
"""
    set_lower_bound(dvar::DecisionVariable, scenario_index::Integer, lower::Number)

Set the lower bound of the scenario-dependent decision variable `dvar` at `scenario_index` to `lower`. If one does not exist, create a new lower bound constraint.
"""
function JuMP.set_lower_bound(dvar::DecisionVariable, scenario_index::Integer, lower::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_lower_bound(dvar, lower)`.")
    new_set = MOI.GreaterThan(convert(Float64, lower))
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if has_lower_bound(dvar, scenario_index)
        # Update existing bound
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dvar).value)
        attr = ScenarioDependentConstraintAttribute(stage(dvar), scenario_index, MOI.ConstraintSet())
        MOI.set(structure(owner_model(dvar)), attr, cindex, new_set)
        # Update proxy
        MOI.set(backend(proxy_), MOI.ConstraintSet(), cindex, new_set)
    else
        # Add new lower bound constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), new_set, stage(dvar), scenario_index)
        # Update proxy if the scenario-dependent decision does not have a lower bound yet
        if !has_lower_bound(DecisionRef(proxy_, index(dvar)))
            MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), new_set)
        end
    end
    return nothing
end
"""
    delete_lower_bound(dvar::DecisionVariable)

Delete the lower bound constraint of the first-stage decision variable `dvar`.
"""
function JuMP.delete_lower_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `delete_lower_bound(dvar, scenario_index)`.")
    JuMP.delete(owner_model(dvar), LowerBoundRef(dvar))
end
"""
    delete_lower_bound(dvar::DecisionVariable, scenario_index::Integer)

Delete the lower bound constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.delete_lower_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `delete_lower_bound(dvar)`.")
    JuMP.delete(owner_model(dvar), LowerBoundRef(dvar), scenario_index)
end
"""
    lower_bound(dvar::DecisionVariable)

Return the lower bound of the first-stage decision variable `dvar`. Error if one does not exist.
"""
function JuMP.lower_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `lower_bound(dvar, scenario_index)`.")
    if !has_lower_bound(dvar)
        error("Decision variable $(dvar) does not have a lower bound.")
    end
    cset = MOI.get(owner_model(dvar), MOI.ConstraintSet(),
                   LowerBoundRef(dvar))::MOI.GreaterThan{Float64}
    return cset.lower
end
"""
    lower_bound(dvar::DecisionVariable, scenario_index::Integer)

Return the lower bound of the scenario-dependent decision variable `dvar` at `scenario_index`. Error if one does not exist.
"""
function JuMP.lower_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `lower_bound(dvar)`.")
    if !has_lower_bound(dvar, scenario_index)
        error("Decision variable $(dvar) at $scenario_index does not have a lower bound.")
    end
    attr = ScenarioDependentConstraintAttribute(stage(dvar), scenario_index, MOI.ConstraintSet())
    cset = MOI.get(owner_model(dvar), attr,
                   LowerBoundRef(dvar))::MOI.GreaterThan{Float64}
    return cset.lower
end
"""
    has_upper_bound(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` has a upper bound.
"""
function JuMP.has_upper_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `upper_bound(dvar, scenario_index)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar))
end
"""
    has_upper_bound(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` has a upper bound at `scenario_index`.
"""
function JuMP.has_upper_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `upper_bound(dvar)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar), scenario_index)
end
"""
    LowerBoundRef(dvar::DecisionVariable)

Return a constraint reference to the upper bound constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.UpperBoundRef(dvar::DecisionVariable)
    moi_ub =  MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}
    return SPConstraintRef{moi_ub, ScalarShape}(owner_model(dvar),
                                                stage(dvar),
                                                moi_ub(index(dvar).value),
                                                ScalarShape())
end
"""
    set_upper_bound(dvar::DecisionVariable, upper::Number)

Set the upper bound of the first-stage decision variable `dvar` to `upper`. If one does not exist, create a new upper bound constraint.
"""
function JuMP.set_upper_bound(dvar::DecisionVariable, upper::Number)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_upper_bound(dvar, scenario_index, upper)`.")
    new_set = MOI.LessThan(convert(Float64, upper))
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if has_upper_bound(dvar)
        # Update existing bound
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dvar).value)
        MOI.set(structure(owner_model(dvar)), MOI.ConstraintSet(), cindex, new_set)
        # Update proxy
        MOI.set(backend(proxy_), MOI.ConstraintSet(), cindex, new_set)
    else
        # Add new upper bound constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), new_set)
        # Update proxy
        MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), new_set)
    end
    return nothing
end
"""
    set_upper_bound(dvar::DecisionVariable, scenario_index::Integer, upper::Number)

Set the upper bound of the scenario-dependent decision variable `dvar` at `scenario_index` to `upper`. If one does not exist, create a new upper bound constraint.
"""
function JuMP.set_upper_bound(dvar::DecisionVariable, scenario_index::Integer, upper::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_upper_bound(dvar, upper)`.")
    new_set = MOI.LessThan(convert(Float64, upper))
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if has_upper_bound(dvar, scenario_index)
        # Update existing bound
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dvar).value)
        attr = ScenarioDependentConstraintAttribute(stage(dvar), scenario_index, MOI.ConstraintSet())
        MOI.set(structure(owner_model(dvar)), attr, cindex, new_set)
        # Update proxy
        MOI.set(backend(proxy_), MOI.ConstraintSet(), cindex, new_set)
    else
        # Add new upper bound constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), new_set, stage(dvar), scenario_index)
        # Update proxy if the scenario-dependent decision does not have a upper bound yet
        if !has_upper_bound(DecisionRef(proxy_, index(dvar)))
            MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), new_set)
        end
    end
    return nothing
end
"""
    delete_upper_bound(dvar::DecisionVariable)

Delete the upper bound constraint of the first-stage decision variable `dvar`.
"""
function JuMP.delete_upper_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `delete_upper_bound(dvar, scenario_index)`.")
    JuMP.delete(owner_model(dvar), UpperBoundRef(dvar))
end
"""
    delete_upper_bound(dvar::DecisionVariable, scenario_index::Integer)

Delete the upper bound constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.delete_upper_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `delete_upper_bound(dvar)`.")
    JuMP.delete(owner_model(dvar), UpperBoundRef(dvar), scenario_index)
end
"""
    upper_bound(dvar::DecisionVariable)

Return the upper bound of the first-stage decision variable `dvar`. Error if one does not exist.
"""
function JuMP.upper_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `upper_bound(dvar, scenario_index)`.")
    if !has_upper_bound(dvar)
        error("Decision $(dvar) does not have a upper bound.")
    end
    cset = MOI.get(owner_model(dvar), MOI.ConstraintSet(),
                   UpperBoundRef(dvar))::MOI.LessThan{Float64}
    return cset.upper
end
"""
    upper_bound(dvar::DecisionVariable, scenario_index::Integer)

Return the upper bound of the scenario-dependent decision variable `dvar` at `scenario_index`. Error if one does not exist.
"""
function JuMP.upper_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `upper_bound(dvar)`.")
    if !has_upper_bound(dvar, scenario_index)
        error("Decision $(dvar) at `scenario_index` does not have a upper bound.")
    end
    attr = ScenarioDependentConstraintAttribute(stage(dvar), scenario_index, MOI.ConstraintSet())
    cset = MOI.get(owner_model(dvar), attr,
                   UpperBoundRef(dvar))::MOI.LessThan{Float64}
    return cset.upper
end
"""
    is_integer(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` is constrained to be integer.
"""
function JuMP.is_integer(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_integer(dvar, scenario_index)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.Integer}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar))
end
"""
    is_integer(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` is constrained to be integer at `scenario_index`.
"""
function JuMP.is_integer(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_integer(dvar)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.Integer}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar), scenario_index)
end
"""
    set_integer(dvar::DecisionVariable)

Add an integrality constraint on the first-stage decision variable `dvar`.
"""
function JuMP.set_integer(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_integer(dvar, scenario_index)`.")
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if is_integer(dvar)
        return nothing
    elseif is_binary(dvar)
        error("Cannot set the decision $(dvar) to integer as it is already binary.")
    else
        # Add new integer constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), MOI.Integer())
        # Update proxy
        MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), MOI.Integer())
    end
    return nothing
end
"""
    set_integer(dvar::DecisionVariable, scenario_index::Integer)

Add an integrality constraint on the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.set_integer(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_integer(dvar)`.")
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if is_integer(dvar, scenario_index)
        return nothing
    elseif is_binary(dvar, scenario_index)
        error("Cannot set the decision $(dvar) at `scenario_index` to integer as it is already binary.")
    else
        # Add new integer constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), MOI.Integer(), stage(dvar), scenario_index)
        # Update proxy if the scenario-dependent decision does not have an integer bound yet
        if !is_integer(DecisionRef(proxy_, index(dvar)))
            MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), MOI.Integer())
        end
    end
    return nothing
end
"""
    unset_integer(dvar::DecisionVariable)

Delete the integrality constraint of the first-stage decision variable `dvar`.
"""
function JuMP.unset_integer(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `unset_integer(dvar, scenario_index)`.")
    JuMP.delete(owner_model(dvar), IntegerRef(dvar))
    return nothing
end
"""
    unset_integer(dvar::DecisionVariable, scenario_index::Integer)

Delete the integrality constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.unset_integer(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `unset_integer(dvar)`.")
    JuMP.delete(owner_model(dvar), IntegerRef(dvar), scenario_index)
    return nothing
end
"""
    IntegerRef(dvar::DecisionVariable)

Return a constraint reference to the integrality constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.IntegerRef(dvar::DecisionVariable)
    moi_int =  MOI.ConstraintIndex{SingleDecision, MOI.Integer}
    return SPConstraintRef{moi_int, ScalarShape}(owner_model(dvar),
                                                 stage(dvar),
                                                 moi_int(index(dvar).value),
                                                 ScalarShape())
end
"""
    is_binary(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` is constrained to be binary.
"""
function JuMP.is_binary(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_binary(dvar, scenario_index)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar))
end
"""
    is_binary(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` is constrained to be binary at `scenario_index`.
"""
function JuMP.is_binary(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_binary(dvar)`.")
    ci = MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), ci, stage(dvar), scenario_index)
end
"""
    set_binary(dvar::DecisionVariable)

Constrain the first-stage decision variable `dvar` to the set ``\\{0,1\\}``.
"""
function JuMP.set_binary(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_binary(dvar, scenario_index)`.")
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if is_binary(dvar)
        return nothing
    elseif is_integer(dvar)
        error("Cannot set the decision $(dvar) to binary as it is already integer.")
    else
        # Add new binary constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), MOI.ZeroOne())
        # Update proxy
        MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), MOI.ZeroOne())
    end
    return nothing
end
"""
    set_binary(dvar::DecisionVariable, scenario_index::Integer)

Constrain the scenario-dependent decision variable `dvar` to the set ``\\{0,1\\}`` at `scenario_index`.
"""
function JuMP.set_binary(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_binary(dvar)`.")
    # Get proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    # Check if bound exists already
    if is_binary(dvar, scenario_index)
        return nothing
    elseif is_integer(dvar, scenario_index)
        error("Cannot set the decision $(dvar) at `scenario_index` to binary as it is already integer.")
    else
        # Add new binary constraint
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), MOI.ZeroOne(), stage(dvar), scenario_index)
        # Update proxy if the scenario-dependent decision does not have a binary bound yet
        if !is_binary(DecisionRef(proxy_, index(dvar)))
            MOI.add_constraint(backend(proxy_), SingleDecision(index(dvar)), MOI.ZeroOne())
        end
    end
end
"""
    unset_binary(dvar::DecisionVariable)

Delete the binary constraint of the first-stage decision variable `dvar`.
"""
function JuMP.unset_binary(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `unset_binary(dvar, scenario_index)`.")
    JuMP.delete(owner_model(dvar), BinaryRef(dvar))
    return nothing
end
"""
    unset_binary(dvar::DecisionVariable, scenario_index::Integer)

Delete the binary constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.unset_binary(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `unset_binary(dvar)`.")
    JuMP.delete(owner_model(dvar), BinaryRef(dvar), scenario_index)
    return nothing
end
"""
    BinaryRef(dvar::DecisionVariable)

Return a constraint reference to the binary constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.BinaryRef(dvar::DecisionVariable)
    moi_bin =  MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}
    return SPConstraintRef{moi_bin, ScalarShape}(owner_model(dvar),
                                                 stage(dvar),
                                                 moi_bin(index(dvar).value),
                                                 ScalarShape())
end
"""
    start_value(dvar::DecisionVariable)

Return the start value of the first-stage decision variable `dvar`.
"""
function JuMP.start_value(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `start_value(dvar, scenario_index)`.")
    return MOI.get(owner_model(dvar), MOI.VariablePrimalStart(), dvar)
end
"""
    start_value(dvar::DecisionVariable, scenario_index::Integer)

Return the start value of the scenario-dependent decision variable `dvar`
at `scenario_index`.
"""
function JuMP.start_value(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `start_value(dvar)`.")
    attr = ScenarioDependentVariableAttribute(stage(dvar), scenario_index, MOI.VariablePrimalStart())
    return MOI.get(owner_model(dvar), attr, dvar)
end
"""
    set_start_value(dvar::DecisionVariable)

Set the start value of the first-stage decision variable `dvar` to `value`.
"""
function JuMP.set_start_value(dvar::DecisionVariable, value::Number)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_start_value(dvar, scenario_index, value)`.")
    MOI.set(backend(proxy(owner_model(dvar), stage(dvar))), MOI.VariablePrimalStart(), index(dvar), Float64(value))
    MOI.set(owner_model(dvar), MOI.VariablePrimalStart(), dvar, Float64(value))
end
"""
    set_start_value(dvar::DecisionVariable, scenario_index::Integer, value::Number)

Set the start value of the scenario-dependent decision variable `dvar`
at `scenario_index` to `value`.
"""
function JuMP.set_start_value(dvar::DecisionVariable, scenario_index::Integer, value::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_start_value(dvar, value)`.")
    attr = ScenarioDependentVariableAttribute(stage(dvar), scenario_index, MOI.VariablePrimalStart())
    MOI.set(owner_model(dvar), attr, dvar, Float64(value))
end

function Base.hash(dvar::DecisionVariable, h::UInt)
    return hash(objectid(owner_model(dvar)), hash(dvar.index, h))
end

JuMP.isequal_canonical(d::DecisionVariable, other::DecisionVariable) = isequal(d, other)
function Base.isequal(dvar::DecisionVariable, other::DecisionVariable)
    return owner_model(dvar) === owner_model(other) && dvar.index == other.index
end

Base.iszero(::DecisionVariable) = false
Base.copy(dvar::DecisionVariable) = DecisionVariable(dvar.stochasticprogram, stage(dvar), dvar.index)
Base.broadcastable(dvar::DecisionVariable) = Ref(dvar)

function DecisionRef(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `DecisionRef(dvar, scenario_index)`.")
    sp = owner_model(dvar)
    return DecisionRef(proxy(sp, stage(dvar)), structure(sp), index(dvar))
end
function DecisionRef(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `DecisionRef(dvar)`.")
    sp = owner_model(dvar)
    return DecisionRef(proxy(sp, stage(dvar)), structure(sp), index(dvar), scenario_index)
end
function DecisionRef(dvar::DecisionVariable, at_stage::Integer, scenario_index::Integer)
    at_stage > stage(dvar) || error("$dvar can only be known after stage $(stage(dvar)).")
    sp = owner_model(dvar)
    return DecisionRef(proxy(sp, at_stage), structure(sp), index(dvar), at_stage, scenario_index)
-end

is_decision_type(::Type{DecisionVariable}) = true

function JuMP.moi_function_type(::Type{DecisionVariable})
    return SingleDecision
end
function JuMP.moi_function_type(::Type{<:Vector{<:DecisionVariable}})
    return VectorOfDecisions
end


text:
using DrWatson
@quickactivate :Swalbe
using CUDA, DataFrames, FileIO
# CUDA.device!(1)

# Fluid dynamics we need for the experiment
"""
    measure_substratewave

Simulates a changing wettability and the response from the film.
"""
function measure_substratewave(
    sys::Swalbe.SysConst,
    device::String;
    h₀ = 1.0,
    ϵ = 0.1,
    wave_x=1,
    wave_y=1,
    sub_speed=100,
    dump = 1000,
    θₛ=ones(sys.Lx, sys.Ly),
    fluid=zeros(sys.Tmax÷dump, sys.Lx*sys.Ly),
    theta=zeros(sys.Tmax÷dump, sys.Lx*sys.Ly),
    dire = "x",
    verbos=true,
    T=Float64
)
    fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py = Swalbe.Sys(sys, device, false, T)
    if device == "CPU"
        for i in 1:sys.Lx, j in 1:sys.Ly
            height[i,j] = h₀ + ϵ * sin(2π*wave_x*(i-1)/sys.Lx) * sin(2π*wave_y*(j-1)/sys.Ly)
        end
    elseif device == "GPU"
        h = zeros(size(height))
        for i in 1:sys.Lx, j in 1:sys.Ly
            h[i,j] = h₀ + ϵ * sin(2π*wave_x*(i-1)/sys.Lx) * sin(2π*wave_y*(j-1)/sys.Ly)
        end
        # theta = CUDA.zeros(Float64, sys.Lx, sys.Ly)
        height = CUDA.adapt(CuArray, h)
    end
    Swalbe.equilibrium!(fout, height, velx, vely, vsq)
    ftemp .= fout
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(height)
            if verbos
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end

        Swalbe.filmpressure!(pressure, height, dgrad, sys.γ, θₛ, sys.n, sys.m, sys.hmin, sys.hcrit)
        Swalbe.∇f!(h∇px, h∇py, pressure, dgrad, height)
        Swalbe.slippage!(slipx, slipy, height, velx, vely, sys.δ, sys.μ)
        # Forces are the pressure gradient and the slippage due to substrate liquid boundary conditions
        Fx .= h∇px .+ slipx
        Fy .= h∇py .+ slipy
        # New equilibrium
        Swalbe.equilibrium!(feq, height, velx, vely, vsq)
        Swalbe.BGKandStream!(fout, feq, ftemp, -Fx, -Fy)
        # New moments
        Swalbe.moments!(height, velx, vely, fout)
        # Measurements, in this case only snapshots of simulational arrays
        Swalbe.snapshot!(fluid, height, t, dumping = dump)
        Swalbe.snapshot!(theta, θₛ, t, dumping = dump)
        move_substrate!(slipx, θₛ, t, sub_speed, direction=dire)
    end
    return fluid, theta
    CUDA.reclaim()
end
# Well not most clean way but one one to freeze the substrate without having the x/0 issue.
function measure_substratewave(
    sys::Swalbe.SysConst,
    device::String,
    move_sub::String;
    h₀ = 1.0,
    ϵ = 0.1,
    wave_x=1,
    wave_y=1,
    sub_speed=100,
    dump=1000,
    θₛ=ones(sys.Lx, sys.Ly),
    fluid=zeros(sys.Tmax÷dump, sys.Lx*sys.Ly),
    theta=zeros(sys.Tmax÷dump, sys.Lx*sys.Ly),
    dire = "x",
    verbos=true,
    T=Float64
)
    println("Simulating a time dependent substrate pattern")
    fout, ftemp, feq, height, velx, vely, vsq, pressure, dgrad, Fx, Fy, slipx, slipy, h∇px, h∇py = Swalbe.Sys(sys, device, false, T)
    if device == "CPU"
        for i in 1:sys.Lx, j in 1:sys.Ly
            height[i,j] = h₀ + ϵ * sin(2π*wave_x*(i-1)/sys.Lx) * sin(2π*wave_y*(j-1)/sys.Ly)
        end
    elseif device == "GPU"
        h = zeros(size(height))
        for i in 1:sys.Lx, j in 1:sys.Ly
            h[i,j] = h₀ + ϵ * sin(2π*wave_x*(i-1)/sys.Lx) * sin(2π*wave_y*(j-1)/sys.Ly)
        end
        # theta = CUDA.zeros(Float64, sys.Lx, sys.Ly)
        height = CUDA.adapt(CuArray, h)
    end
    Swalbe.equilibrium!(fout, height, velx, vely, vsq)
    ftemp .= fout
    for t in 1:sys.Tmax
        if t % sys.tdump == 0
            mass = 0.0
            mass = sum(height)
            if verbos
                println("Time step $t mass is $(round(mass, digits=3))")
            end
        end

        Swalbe.filmpressure!(pressure, height, dgrad, sys.γ, θₛ, sys.n, sys.m, sys.hmin, sys.hcrit)
        Swalbe.∇f!(h∇px, h∇py, pressure, dgrad, height)
        Swalbe.slippage!(slipx, slipy, height, velx, vely, sys.δ, sys.μ)
        # Forces are the pressure gradient and the slippage due to substrate liquid boundary conditions
        Fx .= h∇px .+ slipx
        Fy .= h∇py .+ slipy
        # New equilibrium
        Swalbe.equilibrium!(feq, height, velx, vely, vsq)
        Swalbe.BGKandStream!(fout, feq, ftemp, -Fx, -Fy)
        # New moments
        Swalbe.moments!(height, velx, vely, fout)
        # Measurements, in this case only snapshots of simulational arrays
        Swalbe.snapshot!(fluid, height, t, dumping = dump)
        Swalbe.snapshot!(theta, θₛ, t, dumping = dump)
        if move_sub == "yes"
            # move_substrate!(slipx, θₛ, t, sub_speed, direction=dire)
        end
    end
    return fluid, theta
    CUDA.reclaim()
end

function move_substrate!(θ, input, t, tmove; direction="diagonal")
    if (t % tmove == 0) & (t > 0)
        if direction == "diagonal"
            circshift!(θ, input, (1,1))
        elseif direction == "x"
            circshift!(θ, input, (1,0))
        elseif direction == "y"
            circshift!(θ, input, (0,1))
        end
        input .= θ
    end

    return nothing
end

function pyramidpattern(Lx, Ly; waves=1, mid=1/9, lower=1/18)
    half = Lx÷2
    quat = Lx÷4
    sixt = Lx÷6
    space = zeros(Lx,Ly)
    theta = zeros(Lx,Ly)
    if waves == 1
        for i in 1:Lx÷2
            space[i,i:(half+1-i)] .= (i-1)
            space[half+1-i,i:(half+1-i)] .= (i-1)
            space[i:(half+1-i),i] .= (i-1)
            space[i:(half+1-i),half+1-i] .= (i-1)
        end
        space[half+1:Lx, 1:half] .= - space[1:half,1:half]
        space[1:half, half+1:Lx] .= - space[1:half,1:half]
        space[half+1:Lx, half+1:Lx] .= space[1:half,1:half]
        theta .= mid .+ lower .* space ./ half

    elseif waves == 2
        for i in 1:Lx÷4
            space[i,i:(quat+1-i)] .= (i-1)
            space[quat+1-i,i:(quat+1-i)] .= (i-1)
            space[i:(quat+1-i),i] .= (i-1)
            space[i:(quat+1-i),quat+1-i] .= (i-1)
        end
        space[quat+1:half, 1:quat] .= - space[1:quat,1:quat]
        space[1:quat, quat+1:half] .= - space[1:quat,1:quat]
        space[quat+1:half, quat+1:half] .= space[1:quat,1:quat]

        space[1:half, half+1:Lx] .= space[1:half,1:half]
        space[half+1:Lx, half+1:Lx] .= space[1:half,1:half]
        space[half+1:Lx, 1:half] .= space[1:half,1:half]
        theta .= mid .+ lower .* space ./ quat

    elseif waves == 3
        thr = Lx÷3
        for i in 1:Lx÷6
            space[i,i:(sixt+1-i)] .= (i-1)
            space[sixt+1-i,i:(sixt+1-i)] .= (i-1)
            space[i:(sixt+1-i),i] .= (i-1)
            space[i:(sixt+1-i),sixt+1-i] .= (i-1)
        end
        space[sixt+1:thr, 1:sixt] .= - space[1:sixt,1:sixt]
        space[1:sixt, sixt+1:thr] .= - space[1:sixt,1:sixt]
        space[sixt+1:thr, sixt+1:thr] .= space[1:sixt,1:sixt]

        space[1:thr, thr+1:2*thr] .= space[1:thr,1:thr]
        space[1:thr, 2*thr+1:3*thr] .= space[1:thr,1:thr]
        space[thr+1:2*thr, 1:thr] .= space[1:thr,1:thr]
        space[thr+1:2*thr, thr+1:2*thr] .= space[1:thr,1:thr]
        space[thr+1:2*thr, 2*thr+1:3*thr] .= space[1:thr,1:thr]
        space[2*thr+1:3*thr, 1:thr] .= space[1:thr,1:thr]
        space[2*thr+1:3*thr, thr+1:2*thr] .= space[1:thr,1:thr]
        space[2*thr+1:3*thr, 2*thr+1:3*thr] .= space[1:thr,1:thr]
        theta .= mid .+ lower .* space ./ sixt
    end

    return theta
end

println("Moving Wettability and possible resonaces")
v_lam1_dia = [0, 4900, 490, 49]
v_lam2_dia = [0, 9802, 980, 98]
v_lam3_dia = [0, 14702, 1470, 147]
# Newer simulations
v4 = [19603, 1960, 196]
v5 = [24504, 2450, 245]
v6 = [29404, 2940, 294]
v7 = [34305, 3430, 343]
v8 = [39206, 3920, 392]
v9 = [44107, 4410, 441]
# To pin down the Rayleigh-Plateu instability
v_lam2_dia_more = [490, 245, 164, 123]
for direction in ["diagonal"] #  "diagonal"
    # Different initial volumes
    for waves in [2]# [4,5,6,7,8,9] # 1 2
        n_vel = length(v_lam2_dia_more)
        # speeds = zeros(Int, 4)
        speeds = zeros(Int, n_vel)
        if waves == 1
            speeds .= v_lam1_dia
        elseif waves == 2
            # speeds .= v_lam2_dia
            speeds .= v_lam2_dia_more
        elseif waves == 3
            speeds .= v_lam3_dia
        elseif waves == 4
            speeds .= v4
        elseif waves == 5
            speeds .= v5
        elseif waves == 6
            speeds .= v6
        elseif waves == 7
            speeds .= v7
        elseif waves == 8
            speeds .= v8
        elseif waves == 9
            speeds .= v9
        end
        for speed in speeds # 1 2 3 # [0]
            pattern = "sine"
            ang = 1/9
            println("Simulating moving substrate wettability with pattern $(pattern) and moving direction $(direction) and speed $(speed)")
            sys = Swalbe.SysConst(Lx=512, Ly=512, γ=0.01, δ=1.0, n=3, m=2, hmin=0.07, Tmax=5000000, tdump=5000)
            #sys = Swalbe.SysConst(Lx=512, Ly=512, γ=0.01, δ=1.0, n=3, m=2, hmin=0.07, Tmax=75000, tdump=500)
            df_fluid = Dict()
            df_sub = Dict()
            θₚ = ones(sys.Lx,sys.Ly)
            # Substrate patterning
            if pattern == "sine"
                for i in 1:sys.Lx, j in 1:sys.Ly
                    θₚ[i,j] = ang + 1/18 * sin(2π*waves*(i-1)/sys.Lx) * sin(2π*waves*(j-1)/sys.Ly)
                end
            elseif pattern == "linear"
                θₚ = pyramidpattern(sys.Lx, sys.Ly, waves=waves)
            end
            # Make a cuarray with the substrate pattern
            θ_in = CUDA.adapt(CuArray, θₚ)
            if speed == 0
                # Actual simulation
                fluid, substrate = measure_substratewave(sys, "GPU", "blub", sub_speed=speed, θₛ=θ_in, dire=direction, dump=sys.tdump)
            else
                fluid, substrate = measure_substratewave(sys, "GPU", sub_speed=speed, θₛ=θ_in, dire=direction, dump=sys.tdump)
            end
            println("Writing measurements to Dict")
            # Filling the dataframes
            for t in 1:sys.Tmax÷sys.tdump
                df_fluid["h_$(t*sys.tdump)"] = fluid[t,:]
                df_sub["theta_$(t*sys.tdump)"] = substrate[t,:]
            end
            println("Saving Dict subdirection $direction subvel $speed and $(pattern) $waves to disk")
            save_ang = Int(round(rad2deg(π*ang)))
            save("data/Moving_wettability/height_direc_$(direction)_sp_$(speed)_$(pattern)_$(waves)_$(save_ang)_tmax_$(sys.Tmax)_v2.jld2", df_fluid)
            # bson("data/Moving_wettability/theta_direc_$(direction)_sp_$(speed)_$(pattern)_$(waves)_tmax_$(sys.Tmax)_v2.bson", df_sub)

            CUDA.reclaim()
        end
    end
end

println("Script done, let's have a look at the data :)")


text:
function foo()
    a = 10
    b = 20
    _c1 = a + b
    _c2 = a - b
    _c3 = a * b
    _c4 = a / b
    _c5 = -(a)
    d = 2.0
    _e1 = a + d
    _e2 = a - d
    _e3 = a * d
    _e4 = a / d
    _f = -3.0
    _g = -(a)
end

function add1(x::Int8, y::Int8)::Int16
    return x + y
end

function add2(x::Int16, y::Int16)::Int32
    return x + y
end

function add3(x::Int32, y::Int32)::Int64
    return x + y
end

function add4(x::Int64, y::Int64)::Int64
    return x + y
end

function add5(x::UInt8, y::UInt8)::UInt16
    return x + y
end

function add6(x::UInt16, y::UInt16)::UInt32
    return x + y
end

function add7(x::UInt32, y::UInt32)::UInt64
    return x + y
end

function add8(x::UInt64, y::UInt64)::UInt64
    return x + y
end

function add9(x::Int8, y::UInt16)::UInt32
    return x + y
end

function sub(x::Int8, y::Int8)::Int8
    return x - y
end

function mul(x::Int8, y::Int8)::Int16
    return x * y
end

function fadd1(x::Int8, y::Float64)::Float64
    return x + y
end

function show()
    @assert(fadd1(convert(Int8, 6), 6.0) == 12)
    @assert(add1(convert(Int8, 127), convert(Int8, 1)) == 128)
    @assert(add2(convert(Int16, 32767), convert(Int16, 1)) == 32768)
    @assert(add3(convert(Int32, 2147483647), convert(Int32, 1)) == 2147483648)
    @assert(add4(9223372036854775807, 1) == 9223372036854775808)
    println("OK")
end

if abspath(PROGRAM_FILE) == @__FILE__
    foo()
    show()
end


text:
################################################################################
#
#  Imports
#
################################################################################

import AbstractAlgebra

################################################################################
#
#  Factorization Methods
#
################################################################################

# Square-free Factorization
function factor_sff(f::AbstractAlgebra.Generic.Poly{AbstractAlgebra.gfelem{Int32}})

    # Is monic?

    R = 1
    C = AbstractAlgebra.gcd(f, f')
    w = div(f, C)
end

# Distinct Degree Factorization
function factor_ddf(f::AbstractAlgebra.Generic.Poly{AbstractAlgebra.gfelem{Int32}})

    # Is square-free?

    i = 1
    k = AbstractAlgebra.base_ring(f)

    S, x = AbstractAlgebra.PolynomialRing(k, "x", cached = false)

    fstar = f

    while AbstractAlgebra.degree(fstar) >= 2*i
        r = x^(AbstractAlgebra.characteristic(k))^i - x
        g = AbstractAlgebra.gcd(fstar, r)
        if g != 1
            S = push!(S, (g, i))
            fstar = AbstractAlgebra.divexact(fstar, g)
        end
    end

    if fstar != 1
        S = push!(S, (fstar, AbstractAlgebra.degree(fstar)))
    elseif S == []
        return [(f, 1)]
    else
        return S
    end
end


text:
using SQLite

name = "Smith, Steve"
jerseys =  Dict("Smith, Steve" => 99)
sqlbool(tf::Bool) = if(tf) "TRUE" else "FALSE" end

db = SQLite.DB() # no filename given, so create an in-memory temporary
SQLite.execute!(db, "create table players (id integer primary key,
                                           name text,
                                           score number,
                                           active bool,
                                           jerseynum integer)")

SQLite.query(db, "INSERT INTO players (name, score, active, jerseynum) values ('Jones, James', 9, 'FALSE', 99)")
SQLite.query(db, "UPDATE players SET name = ?, score = ?, active = ? WHERE jerseynum = ?";
    values = ["Smith, Steve", 42, sqlbool(true), jerseys[name]])

tbl = SQLite.query(db, "SELECT * from players")
println(tbl)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export BVector

function make_test_bvector(name:: Symbol;
    isCalibrated :: Bool = true,
    increasing :: Symbol = :increasing, offset = 0.0)

    n = 4;
    defaultValue = collect(LinRange(0.3, 0.2, n));
    ca = BVector(name, string(name), string(name), defaultValue .+ 0.1, defaultValue,
    -2.0 + offset, 3.0 + offset, isCalibrated,
    increasing);
    return ca
end

## -------------  Basics

# Need those because they are otherwise defined in terms of `default_value`
Base.size(ca :: BVector{T1}) where T1 = size(ca.defaultDxV);
Base.length(ca :: BVector{T1}) where T1 = length(ca.defaultDxV);
Base.show(io :: IO,  p :: BVector{T1}) where T1 =
    print(io, "BVector:  " * show_string(p));

# is_calibrated(ca :: BVector{T1}) where T1 = ca.isCalibrated;

calibrated_value(ca :: BVector{T1}) where T1 =
    is_calibrated(ca)  ?  ca.dxV  :  nothing;
calibrated_lb(ca :: BVector{T1}) where T1 = zeros(T1, length(ca));
calibrated_ub(ca :: BVector{T1}) where T1 = ones(T1, length(ca));

function validate(ca :: BVector{T1}; silent :: Bool = false) where T1
    isValid = all(x -> x >= zero(T1), ca.dxV)  &&  all(x -> x <= one(T1), ca.dxV);
    return isValid
end

is_increasing(iv :: BVector{T1}) where T1 = (iv.increasing == :increasing);
is_decreasing(iv :: BVector{T1}) where T1 = (iv.increasing == :decreasing);
is_nonmonotone(iv :: BVector{T1}) where T1 = (iv.increasing == :nonmonotone);
# For consistency, this returns a Vector
lb(iv :: BVector{T1}) where T1 = fill(iv.lb, size(iv));
ub(iv :: BVector{T1}) where T1 = fill(iv.ub, size(iv));
# Scalar bounds across all elements
scalar_lb(iv :: BVector{T1}) where T1 = iv.lb;
scalar_ub(iv :: BVector{T1}) where T1 = iv.ub;

## -----------  Retrieve

"""
	$(SIGNATURES)

Returns all values of a `BoundedVector`.
"""
value(iv :: BVector{T1}) where T1 =
    dx_to_values(iv, iv.dxV);

default_value(iv :: BVector{T1}) where T1 =
    dx_to_values(iv, iv.defaultDxV);

# Input is in terms of untransformed units.
function fix!(iv :: BVector{T1}; pValue = nothing) where T1
   iv.isCalibrated = false;
   isnothing(pValue)  ||  set_default_value!(iv, pValue);
end

# The input is in untransformed units.
function set_value!(iv :: BVector{T1}, vIn;
    skipInvalidSize = false) where T1

    oldValue = value(iv);
    if size(iv) == size(vIn)
        iv.dxV = values_to_dx(iv, vIn);
    else
        @warn("""
            Wrong size for $iv
            Given: $(size(vIn))
            Expected: $(size(iv))
            """);
        if !skipInvalidSize
            error("Stopped");
        end
    end
    return oldValue
end

# Input is in terms of dxV
function set_calibrated_value!(iv :: BVector{T1}, vIn :: AbstractVector{T1};
    skipInvalidSize = false) where T1

    if size(iv) == size(vIn)
        iv.dxV .= vIn;
    else
        @warn """
            Invalid size of new calibrated values for $iv
            Given: $(size(vIn))
            Expected: $(size(iv))
            """
        skipInvalidSize  ||  error("Stopped");
    end
end

# Input is in untransformed units.
function set_default_value!(iv :: BVector{T1}, vIn :: AbstractVector{T1}) where T1
    iv.defaultDxV = values_to_dx(iv, vIn);
end

# Bounds have to be scalar or vector with all equal elements.
function set_bounds!(iv :: BVector{T1}; lb = nothing, ub = nothing)  where T1
    isnothing(lb)  ||  set_lower_bound!(iv, lb);
    isnothing(ub)  ||  set_upper_bound!(iv, ub);
    @assert size(ModelParams.lb(iv)) == size(ModelParams.ub(iv)) == size(iv);
end

set_lower_bound!(iv :: BVector{T1}, lb) where T1 = iv.lb = vec_to_scalar(lb);
set_upper_bound!(iv :: BVector{T1}, ub) where T1 = iv.ub = vec_to_scalar(ub);

vec_to_scalar(v :: F1) where F1 <: Real = v;

function vec_to_scalar(v :: AbstractVector{F1}) where F1 <: Real
    @assert all(isapprox.(v, first(v)))  "All elements expected to be the same.";
    return first(v)
end

# -------------


text:
using PrettyEinsum
using Test

using TensorOperations

@testset "PrettyEinsum.jl" begin
    rsl1 = zeros(ComplexF64, 5, 5, 2)
    rsl2 = zeros(ComplexF64, 5, 5, 2)
    rsl4 = zeros(ComplexF64, 5, 5, 2)

    l = randn(ComplexF64, 5, 5)
    a = randn(ComplexF64, 5, 5, 2)
    r = randn(ComplexF64, 5, 5)
    H = randn(ComplexF64, 2, 2, 2, 2)

    @tensoropt !(p, p_lower_r, p_upper_l, p_upper_r) rsl1[l, r, p] =
        l[l, L] *
        a[L, CL, p_upper_l] *
        a[CL, CR, p_upper_r] *
        r[CR, R] *
        conj(a[r, R, p_lower_r]) *
        H[p, p_lower_r, p_upper_l, p_upper_r]

    t, _ = @optimalcontractiontree !(p, p_lower_r, p_upper_l, p_upper_r) rsl1[l, r, p] =
        l[l, L] *
        a[L, CL, p_upper_l] *
        a[CL, CR, p_upper_r] *
        r[CR, R] *
        conj(a[r, R, p_lower_r]) *
        H[p, p_lower_r, p_upper_l, p_upper_r]
    @test t == [[1, [6, 2]], [5, [3, 4]]]

    # inplace
    Σein!"""
    diagram:
    #      ______[2|l,r,p]____(r@2=l@3)____[3|l,r,p]_
    #     |        |                         |       |
    # (r@1=l@2)  (p@2=p_up_L@6)   (p@3=p_up_R@6) (r@3=l@4)
    #   [1|l,r]  [6|p_lo_L,p_lo_R,p_up_L,p_up_R]   [4|l,r]
    #    l@1     p_lo_L@6         (p@5=p_lo_R@6) (r@4=r@5)
    #                                        |       |
    #                               l@5___[5*|l,r,p]_| -> [>|l@1,l@5,p_lo_L@6]
    contraction: [(1,(6,2)),(5,(3,4))]
    """(rsl2, l, a, a, r, a, H)

    # allocate
    rsl3 = Σein"""
    diagram:
    #      ______[2|l,r,p]____(r@2=l@3)____[3|l,r,p]_
    #     |        |                         |       |
    # (r@1=l@2)  (p@2=p_up_L@6)   (p@3=p_up_R@6) (r@3=l@4)
    #   [1|l,r]  [6|p_lo_L,p_lo_R,p_up_L,p_up_R]   [4|l,r]
    #    l@1     p_lo_L@6         (p@5=p_lo_R@6) (r@4=r@5)
    #                                        |       |
    #                               l@5___[5*|l,r,p]_| -> [>|l@1,l@5,p_lo_L@6]
    contraction: [(1,(6,2)),(5,(3,4))]
    """(l, a, a, r, a, H)

    foo!(x, A, L, R, h) = @Σein! """
    diagram:
    #       __[$A@2|l,r,p]____(r@2=l@3)____[$A@3|l,r,p]__
    #      |       |                            |        |
    #  (r@1=l@2) (p@2=p_up_L@6)      (p@3=p_up_R@6)  (r@3=l@4)
    # [$L@1|l,r] [$h@6|p_lo_L,p_lo_R,p_up_L,p_up_R] [$R@4|l,r]
    #     l@1    p_lo_L@6            (p@5=p_lo_R@6)  (r@4=r@5)
    #                                           |        |
    #                               l@5___[$A@5*|l,r,p]__| -> [$x@>|l@1,l@5,p_lo_L@6]
    contraction: [(1,(6,2)),(5,(3,4))]
    """
    foo!(rsl4, a, l, r, H)

    foo(A, L, R, h) = @Σein """
    diagram:
    #       __[$A@2|l,r,p]____(r@2=l@3)____[$A@3|l,r,p]__
    #      |       |                            |        |
    #  (r@1=l@2) (p@2=p_up_L@6)      (p@3=p_up_R@6)  (r@3=l@4)
    # [$L@1|l,r] [$h@6|p_lo_L,p_lo_R,p_up_L,p_up_R] [$R@4|l,r]
    #     l@1    p_lo_L@6            (p@5=p_lo_R@6)  (r@4=r@5)
    #                                           |        |
    #                               l@5___[$A@5*|l,r,p]__| -> [>|l@1,l@5,p_lo_L@6]
    contraction: [(1,(6,2)),(5,(3,4))]
    """
    rsl5 = foo(a, l, r, H)

    # @test rsl1 == rsl2 == rsl3 == rsl4 = rsl5
    @test rsl1 == rsl2
    @test rsl1 == rsl3
    @test rsl1 == rsl4
    @test rsl1 == rsl5
end