text:
using RobustOptimization
using JuMP
using Test
using MathOptInterface
const MOI = MathOptInterface
@testset "RobustOptimization.jl" begin
    # Write your tests here.
    @testset "RobustOptimization.jl" begin
    # Write your tests here.

    m = RobustModel()

    @variable(m, 0 <= y <= 1,Uncertain(1))

    @variable(m, 1 <= x_interval <= 3.5,Decision(1,"test",m.uncertainVariables))
    #Add VariableTest
    @testset "variables" begin
        #init
        mod = RobustModel()
        #add vars
        @variable(mod, 0 <= y <= 1,Uncertain(1))
        @variable(mod, 0 <= z <= 1,Uncertain(1))
        @variable(mod, 1 <= x_interval <= 3.5,Decision(1,"test",mod.uncertainVariables))

        #test add vars
        @test mod.var_to_name[3] == "x_interval"
        @test mod.var_to_name[1] == "y"

        #test infos
        @test has_lower_bound(x_interval) == true
        @test lower_bound(x_interval) == 1.0
        @test has_upper_bound(x_interval) == true
        @test upper_bound(x_interval)== 3.5

        #test setter
        set_lower_bound(x_interval,0.5)
        @test lower_bound(x_interval) == 0.5
        set_upper_bound(x_interval,5.5)
        @test upper_bound(x_interval) == 5.5

        #test start value
        set_start_value(x_interval,2.0)
        @test start_value(x_interval) == 2.0

        #test delete
        delete_lower_bound(x_interval)
        @test has_lower_bound(x_interval) == false
        delete_upper_bound(x_interval)
        @test has_upper_bound(x_interval) == false

        #test fixing
        fix(x_interval, 4)
        @test is_fixed(x_interval) == true

        #test unfix
        unfix(x_interval)
        @test is_fixed(x_interval) == false

        #test binary
        set_binary(x_interval)
        @test is_binary(x_interval) == true
        unset_binary(x_interval)
        @test is_binary(x_interval) == false

        #test integer
        set_integer(x_interval)
        @test is_integer(x_interval) == true
        unset_integer(x_interval)
        @test is_integer(x_interval) == false

        #test delete
        decLen = length(mod.decisionVariables)
        delete(mod,x_interval)
        @test length(mod.decisionVariables)-decLen == -1

    end

    @testset "constraints" begin
        model = RobustModel()
        @variable(model, 0 <= y <= 1,Uncertain(1))
        @constraint(model, uncConst, 1.5y >= 1, UncertaintySetConstraint("1"))
        @constraint(model, modConst, 2y <= 2, UncertainConstraint("2"))
        @constraint(model, modConstt, 2y <= 2, UncertainConstraint("2"))
        @test name(uncConst) == "uncConst"
        print(model,"\n")

        @test string(model.uncertainConstraints[uncConst.type_idx].func)== "2 y - 2"
        @test typeof(model.uncertainConstraints[uncConst.type_idx].set) == MathOptInterface.LessThan{Float64}
        #JuMP.coefficient() #model.uncertainConstraints[uncConst.type_idx], model.uncertainVariables[y.type_idx], 4)
        print(typeof(model.uncertainConstraints[uncConst.type_idx].func),"\n")

    end

    #@constraint(m, c1, 2*x_interval + 2*y_interval <= 10)
    #@objective(m, Max, x_interval + y_interval)
    @test m.var_to_name[2] == "x_interval"
    #TODO get inner of Decision
    @test m.var_to_name[1] == "y"
    @test m.nextvaridx == length(m.uncertainVariables)+length(m.decisionVariables)
    #con =m.uncertainConstraints[RobustOptimization.Int(1)]
    #@test con.set == MathOptInterface.LessThan{Float64}(10.0)

    ########### TEST Variable Basics

    print(m.objectivesense)
end

end


text:
using Plots
using ForwardDiff
using LinearAlgebra

a = 1.0
b = 100.0

f(x) = (a - x[1])^2.0 + b * (x[2] - x[1]^2.0)^2.0
f(x, y) = f([x; y])
∇f(x) = ForwardDiff.gradient(f, x)
∇²f(x) = ForwardDiff.hessian(f, x)

n = 1000
x_limits = range(-1.0, stop = 1.0, length = n)
y_limits = range(-1.0, stop = 1.0, length = n)

x = zeros(2)

function solve(x;
        max_iter = 100,
        visualize = true)
    visualize && (plt = Plots.contour(x_limits, y_limits, f))
    visualize && (plt = plot!((x[1], x[2]), markershape = :circle, color = :cyan, outline = :cyan))

    for i = 1:max_iter
        grad = ∇f(x)
        println("x:  $x")
        println("grad norm: $(norm(grad))")
        norm(grad) < 1.0e-3 && break

        hess = ∇²f(x)
        # Δ = hess \ grad
        Δ = grad

        α = 1.0
        iter = 1
        while norm(∇f(x - α * Δ)) > norm(grad) && iter < 100
            α *= 0.5
            iter += 1
        end

        x .-= α * Δ

        if i % 10 == 0 && visualize
            plt = plot!((x[1], x[2]),
                markershape = :circle, color = :cyan, outline = :cyan,
                legend = false)
            display(plt)
        end
    end
    print("iters: $i")
    return x
end

solve([0.0, -1.0], max_iter = 1000)

ext: jl

max_issues_repo_name: thowell/DirectMotionPlanning


text:
module NestedArrays

import Base: @propagate_inbounds, @_inline_meta
import Base: IndexStyle, size, indices, getindex, setindex!, transpose, ctranspose
import Base.Broadcast: _containertype, promote_containertype, broadcast_c, _broadcast!

export NestedArray, NestedVector, NestedMatrix, nested

"""
    NestedArray(a)

Construct a `NestedArray` wrapping array `a`. The nested array applies broadcasting and
linear algebra operations in a nested sense. The elements of `a` are expected to behave
similarly to `AbstractArray`, but no type requirement is imposed.
"""
struct NestedArray{T, n, A <: AbstractArray{T, n}} <: AbstractArray{T, n}
    a::A
end

const NestedVector{T, A <: AbstractVector{T}} = NestedArray{T, 1, A}
const NestedMatrix{T, A <: AbstractVector{T}} = NestedArray{T, 2, A}
const NestedVecOrMat{T, A <: AbstractVecOrMat{T}} = Union{NestedVector{T, A}}

NestedArray(a::AbstractArray{T, n}) where {T, n} = NestedArray{T, n, typeof(a)}(a)
NestedArray{T}(a::AbstractArray{T, n}) where {T, n} = NestedArray{T, n, typeof(a)}(a)
NestedArray{T, n}(a::AbstractArray{T, n}) where {T, n} = NestedArray{T, n, typeof(a)}(a)
NestedVector(a::AbstractVector{T}) where {T} = NestedArray{T, 1, typeof(a)}(a)
NestedMatrix(a::AbstractMatrix{T}) where {T} = NestedArray{T, 2, typeof(a)}(a)

"""
    nested(a::AbstractArray)

Construct a `NestedArray` wrapping `a`, which allows for recursive broadcasting, transpose
and linear algebra.
"""
nested(a::AbstractArray) = NestedArray(a)

@inline size(a::NestedArray) = size(a.a)
@inline indices(a::NestedArray) = indices(a.a)
@inline IndexStyle(a::NestedArray) = IndexStyle(a.a)
@propagate_inbounds getindex(a::NestedArray, i...) = getindex(a.a, i...)
@propagate_inbounds setindex!(a::NestedArray, v, i...) = setindex!(a.a, v, i...)

# Transposition is recursive (same as Base, at the moment...)
# TODO - make (c)transpose(::AbstractMatrix) lazy
#      - make (c)transpose(::AbstractVecOrMat) non-recursive
#      - make (c)transpose(::NestedMatOrVec) recursive
transpose(a::NestedVecOrMat) = NestedArray(transpose(a.a))
ctranspose(a::NestedVecOrMat) = NestedArray(ctranspose(a.a))

# Broadcast results in a NestedArray and is recursive
_containertype(::Type{<:NestedArray}) = NestedArray
promote_containertype(::Type{NestedArray}, ::Type{Any}) = NestedArray
promote_containertype(::Type{Any}, ::Type{NestedArray}) = NestedArray
promote_containertype(::Type{NestedArray}, ::Type{Array}) = NestedArray
promote_containertype(::Type{Array}, ::Type{NestedArray}) = NestedArray
# TODO: tuples and sparse arrays

@inline function broadcast_c(f, ::Type{NestedArray}, a::NestedArray)
    NestedArray(broadcast(x -> (@_inline_meta; broadcast(f, x)), a.a))
end

@inline function broadcast_c(f, ::Type{NestedArray}, a1::NestedArray, a2::NestedArray)
    NestedArray(broadcast((x1, x2) -> (@_inline_meta; broadcast(f, x1, x2)), a1.a, a2.a))
end

# TODO broken for a1::Array?
@inline function broadcast_c(f, ::Type{NestedArray}, a1, a2::NestedArray)
    NestedArray(broadcast((x1, x2) -> (@_inline_meta; broadcast(y2 -> (@_inline_meta; f(x1, y2)), x2)), a1, a2.a))
end

# TODO broken for a2::Array?
@inline function broadcast_c(f, ::Type{NestedArray}, a1::NestedArray, a2)
    NestedArray(broadcast((x1, x2) -> (@_inline_meta; broadcast(y1 -> (@_inline_meta; f(y1, x2)), x1)), a1.a, a2))
end

# TODO finish making this generic
# TODO seperate `all(nest)` branch using dispatch instead?
@generated function broadcast_c(f, ::Type{NestedArray}, a...)
    nest = map(T -> T <: NestedArray, a.params)
    n = length(nest)
    if all(nest)
        expr_as = [:(a[$i]) for i = 1:n]
        expr_xs = [Symbol("x$i") for i = 1:n]
        return quote
            @_inline_meta
            NestedArray(broadcast(f, tuple($(expr_xs...)) -> broadcast(f, $(expr_xs...)), $(expr_as...)))
        end
    else
        error("Generic mixed nested/non-nested broadcasting is not implemented yet")
    end
end

# Same for broadcast!
# TODO: eliminating inner allocations
@inline function broadcast_c!(f, ::Type{NestedArray}, c, a::NestedArray)
    broadcast!(x -> (@_inline_meta; broadcast(f, x)), c, a.a)
    return c
end

@inline function broadcast_c!(f, ::Type{NestedArray}, c, a1::NestedArray, a2::NestedArray)
    NestedArray(broadcast!((x1, x2) -> (@_inline_meta; broadcast(f, x1, x2)), c, a1.a, a2.a))
    return c
end

# TODO broken for a1::Array?
@inline function broadcast_c!(f, ::Type{NestedArray}, c, a1, a2::NestedArray)
    NestedArray(broadcast!((x1, x2) -> (@_inline_meta; broadcast(y2 -> (@_inline_meta; f(x1, y2)), x2)), c, a1, a2.a))
    return c
end

# TODO broken for a2::Array?
@inline function broadcast_c!(f, ::Type{NestedArray}, c, a1::NestedArray, a2)
    NestedArray(broadcast!((x1, x2) -> (@_inline_meta; broadcast(y1 -> (@_inline_meta; f(y1, x2)), x1)), c, a1.a, a2))
    return c
end

# TODO finish making this generic
@generated function broadcast_c!(f, ::Type{NestedArray}, c, a...)
    nest = map(T -> T <: NestedArray, a.params)
    n = length(nest)
    if all(nest)
        expr_as = [:(a[$i].a) for i = 1:n]
        expr_xs = [Symbol("x$i") for i = 1:n]
        return quote
            @_inline_meta
            broadcast!(f, tuple($(expr_xs...)) -> broadcast(f, $(expr_xs...)), c, $(expr_as...))
            return c
        end
    else
        error("Generic mixed nested/non-nested broadcasting is not implemented yet")
    end
end

end # module


text:
"""
Test an idea about a variant of iterated gap trees:
Include the non-tree edges by rotate the capacities in the associated cycles.
"""
module CycleGap

import GraphIdx
import GraphIdx: Graph, EdgeGraph, Edge, Weights, WeightedGraph
import GraphIdx: PrimMstMem, prim_mst_edges
import GraphIdx.Tree: RootedTree

import TreeLas.TreeDP: TreeDPMem, tree_dp!
import TreeLas.Dual: dual!
import TreeLas.MGT: update_tree!
import TreeLas: Cycle
import TreeLas.Cycle: CycleBasis
import TreeLas: GapLas
import TreeLas.GapLas: GapMem, gaplas, Sol

"""
Memory to compute an iteration using the cycle shift idea.

!!! note

    At the moment, CyclceBasis does not have a `reset!` method or alike, so this has to be
    `mutable struct`

"""
mutable struct CycMem{N, WL<:Weights{Float64}}
    sol::Sol{N}
    gamma::Vector{Float64}
    gap_mem::GapMem{N, WL}
    cycles::CycleBasis
    λ_tree::Vector{Float64}
end

function CycMem(y::Array{Float64,N}, graph::Graph, lambda::W) where {N, W<:Weights{Float64}}
    gap_mem = GapMem(y, graph, lambda)
    λ_tree = collect(lambda, GraphIdx.num_nodes(graph))
    return CycMem(
        gap_mem.sol,
        gap_mem.gamma,
        gap_mem,
        CycleBasis(GraphIdx.num_edges(graph), GraphIdx.num_nodes(graph)),
        λ_tree,
    )
end

function cyclas(
    y::Array{Float64, N},
    graph::Graph,
    lambda::Weights{Float64},
    mu::Weights{Float64} = GraphIdx.Ones{Float64}();
    max_iter::Int = 5,
    verbose::Bool = true,
) where {N}
    gaplas(y, graph, lambda, mu, CycMem, max_iter = max_iter, verbose = verbose)
end

function GapLas.gaplas!(
    cmem::CycMem{N, W1},
    y::Array{Float64,N},
    graph::GraphT,
    lambda::Weights{Float64},
    mu::W2,
) where {N, GraphT<:Graph, W1<:Weights{Float64}, W2<:Weights{Float64}}
    mem = cmem.gap_mem
    GapLas.find_gap_tree!(mem, y, graph, lambda)

    cmem.cycles = Cycle.CycleBasis(graph, mem.mst.parent)
    tlam0, tlam = Cycle.extract_rotate(cmem.cycles, lambda)
    cmem.λ_tree .= tlam
    length(tlam) < 30 && @debug tlam

    tree_dp!(mem.sol.x, mem.y_tree, mem.tree, GraphIdx.Vec(cmem.λ_tree), mu, mem.dp_mem)
    mem.α_tree .= vec(mem.sol.x) .- vec(mem.y_tree)
    dual!(mem.α_tree, mem.dp_mem.proc_order, mem.mst.parent)

    fix_alpha!(mem.α_tree, tlam0, mem.sol.x, mem.mst.parent)

    update_tree!(
        mem.sol.α,
        mem.α_tree,
        mem.mst.selected,
        mem.egraph,
        mem.mst.parent,
    )
end

function fix_alpha!(alpha, lam, x, pi)
    n = length(alpha)
    n < 30 && @debug "before" alpha lam x
    for i = 1:n
        if i == pi[i]
            continue
        end
        let ai = alpha[i]
            if abs(ai) > lam[i]
                alpha[i] = sign(ai) * lam[i]
                let diff = alpha[i] - ai
                    x[i] += diff
                    x[pi[i]] -= diff
                end
            end
        end
    end
    n < 30 && @debug "after" alpha lam x
end

end # CycleGap

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
"""
environment.jl

This file defines the structs that capture information about
the experiment such as `FilteringEnvironment`,
`MergingEnvironment` and `MergingEnvironmentLower`
"""

# The upper merge
"""
    MergingEnvironment

Road network with merge point. This is for the upper merge i.e. a and b1

# Internals

- `roadway::Roadway{Float64}` contains all the road segment and lane information
- `merge_point::VecSE2{Float64}` coordinate of the merge point in cartesian frame (0.0, 0.0, 0.0) by default
- `merge_proj::RoadProjection{Int64, Float64}` projection of the merge point on the roadway
- `merge_index::RoadIndex`

# Example: Caution
- If you use all default values while constructing, does not work, so I provide merge_point explicitly
```julia
env_interaction = MergingEnvironment(merge_point = VecSE2(1064.5227,959.1559,-2.8938))
```
"""
@with_kw struct MergingEnvironment
    roadway::Roadway{Float64} = make_roadway_interaction_with_extensions()
    merge_point::VecSE2{Float64} = VecSE2(1064.5227, 959.1559, -2.8938)
    merge_proj::RoadProjection{Int64, Float64} = proj(merge_point, roadway)
    merge_index::RoadIndex{Int64, Float64} = RoadIndex(merge_proj.curveproj.ind, merge_proj.tag)
end

"""
    main_lane(env::MergingEnvironment)
returns the main lane of the merging scenario. Hardcoded to `b1`
"""
main_lane(env::MergingEnvironment) = env.roadway[LaneTag(1, 2)]
main_lane_tag(env::MergingEnvironment) = LaneTag(1,2)

"""
    merge_lane(env::MergingEnvironment)
returns the merging lane of the merging scenario. Harcoded to `a`
"""
merge_lane(env::MergingEnvironment) = env.roadway[LaneTag(1, 1)]
merging_lane_tag(env::MergingEnvironment) = LaneTag(1,1)

# Lower merge environment
"""
struct MergingEnvironmentLower
Environment for the lower merge i.e. road segment 4, lane tag 1 and lane tag 2

# Example
```julia
c = CooperativeIDM(env=MergingEnviornmentLower())
```
"""
@with_kw struct MergingEnvironmentLower
    roadway::Roadway{Float64} = make_roadway_interaction_with_extensions()
    merge_point::VecSE2{Float64} = VecSE2(1110.827, 944.598, 0.214)
    merge_proj::RoadProjection{Int64, Float64} = proj(merge_point, roadway)
    merge_index::RoadIndex{Int64, Float64} = RoadIndex(merge_proj.curveproj.ind, merge_proj.tag)
end

"""
    main_lane(env::MergingEnvironmentLower)
returns the main lane of the merging scenario. Hardcoded to `f1`
"""
main_lane(env::MergingEnvironmentLower) = env.roadway[LaneTag(4, 2)]

"""
main_lane_tag(env::MergingEnvironmentLower) = LaneTag(4,2)

Returns the lane tag of the main lane i.e. f1
"""
main_lane_tag(env::MergingEnvironmentLower) = LaneTag(4,2)

"""
    merge_lane(env::MergingEnvironmentLower)
returns the merging lane of the merging scenario. Harcoded to `g`
"""
merge_lane(env::MergingEnvironmentLower) = env.roadway[LaneTag(4, 1)]

"""
merge_lane_tag(env::MergingEnvironmentLower) = LaneTag(4,1)

Returns the lane tag of the merging lane. Used by the `find_merge_vehicle` method
in cooperative IDM
"""
merging_lane_tag(env::MergingEnvironmentLower) = LaneTag(4,1)

# struct: FilteringEnvironment
"""
@with_kw struct FilteringEnvironment

- Contains the roadway and trajectory information
- To enable not having to pass roadway and traj as input args to all functions

# Fields
- roadway
- traj
- timestep
- `mergeenv`: Upper merging section (`MergingEnvironment`) vs lower (`MergingEnvironmentLower`)

# Example
```julia
f = FilteringEnvironment()
```
"""
@with_kw struct FilteringEnvironment
    roadway::Roadway{Float64} = make_roadway_interaction_with_extensions()
    traj = read_veh_tracks(roadway=roadway)
    timestep::Float64 = 0.1
    mergeenv = MergingEnvironment()
end

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
@testset "Test Symbol Table" begin
    @testset "Scope" begin
        expected = Dict(
            "a" => m.MonkeySymbol("a", m.GlobalScope, 0, nothing),
            "b" => m.MonkeySymbol("b", m.GlobalScope, 1, nothing),
            "c" => m.MonkeySymbol("c", m.LocalScope, 0, nothing),
            "d" => m.MonkeySymbol("d", m.LocalScope, 1, nothing),
            "e" => m.MonkeySymbol("e", m.LocalScope, 0, nothing),
            "f" => m.MonkeySymbol("f", m.LocalScope, 1, nothing),
        )

        g = m.SymbolTable()

        a = m.define!(g, "a")
        @test a == expected["a"]

        b = m.define!(g, "b")
        @test b == expected["b"]

        l1 = m.SymbolTable(; outer = g)

        c = m.define!(l1, "c")
        @test c == expected["c"]

        d = m.define!(l1, "d")
        @test d == expected["d"]

        l2 = m.SymbolTable(; outer = l1)

        e = m.define!(l2, "e")
        @test e == expected["e"]

        f = m.define!(l2, "f")
        @test f == expected["f"]

        @test m.resolve(g, "a")[1] == expected["a"]
        @test m.resolve(g, "b")[1] == expected["b"]

        @test m.resolve(l1, "a")[1] == expected["a"]
        @test m.resolve(l1, "b")[1] == expected["b"]
        @test m.resolve(l1, "c")[1] == expected["c"]
        @test m.resolve(l1, "d")[1] == expected["d"]

        @test m.resolve(l2, "a")[1] == expected["a"]
        @test m.resolve(l2, "b")[1] == expected["b"]
        @test m.resolve(l2, "e")[1] == expected["e"]
        @test m.resolve(l2, "f")[1] == expected["f"]
    end

    @testset "Shadowing Variables" begin
        g = m.SymbolTable()
        a = m.define!(g, "a")
        b = m.define!(g, "b")

        l1 = m.SymbolTable(; outer = g)
        @test m.resolve(l1, "a")[1] == m.MonkeySymbol("a", m.GlobalScope, 0, nothing)

        a1 = m.define!(l1, "a")
        @test a1 == m.MonkeySymbol("a", m.LocalScope, 0, nothing)
        @test m.resolve(l1, "a")[1] == m.MonkeySymbol("a", m.LocalScope, 0, nothing)
        @test m.resolve(l1, "b")[1] == m.MonkeySymbol("b", m.GlobalScope, 1, nothing)

        l2 = m.SymbolTable(; outer = l1, within_loop = true)
        @test m.resolve(l2, "a")[1] ==
              m.MonkeySymbol("a", m.OuterScope, 0, m.SymbolPointer(1, m.LocalScope, 0))
        @test m.resolve(l2, "b")[1] == m.MonkeySymbol("b", m.GlobalScope, 1, nothing)

        l3 = m.SymbolTable(; outer = l2, within_loop = true)
        @test m.resolve(l3, "a")[1] ==
              m.MonkeySymbol("a", m.OuterScope, 0, m.SymbolPointer(2, m.LocalScope, 0))
        @test m.resolve(l3, "b")[1] == m.MonkeySymbol("b", m.GlobalScope, 1, nothing)
    end

    @testset "Builtins" begin
        expected = Dict(
            "a" => m.MonkeySymbol("a", m.BuiltinScope, 0, nothing),
            "c" => m.MonkeySymbol("c", m.BuiltinScope, 1, nothing),
            "e" => m.MonkeySymbol("e", m.BuiltinScope, 2, nothing),
            "f" => m.MonkeySymbol("f", m.BuiltinScope, 3, nothing),
        )

        g = m.SymbolTable()
        l1 = m.SymbolTable(; outer = g)
        l2 = m.SymbolTable(; outer = l1)

        for (i, k) in expected |> keys |> collect |> sort |> enumerate
            m.define_builtin!(g, k, i - 1)
        end

        for s in [g, l1, l2]
            for (k, v) in expected
                @test m.resolve(s, k)[1] == v
            end
        end
    end

    @testset "Free Variables" begin
        g = m.SymbolTable()
        m.define!(g, "a")
        m.define!(g, "b")

        l1 = m.SymbolTable(; outer = g)
        m.define!(l1, "c")
        m.define!(l1, "d")

        l2 = m.SymbolTable(; outer = l1)
        m.define!(l2, "e")
        m.define!(l2, "f")

        for (table, expected_symbols, expected_free_symbols) in [
            (
                l1,
                [
                    m.MonkeySymbol("a", m.GlobalScope, 0, nothing),
                    m.MonkeySymbol("b", m.GlobalScope, 1, nothing),
                    m.MonkeySymbol("c", m.LocalScope, 0, nothing),
                    m.MonkeySymbol("d", m.LocalScope, 1, nothing),
                ],
                [],
            ),
            (
                l2,
                [
                    m.MonkeySymbol("a", m.GlobalScope, 0, nothing),
                    m.MonkeySymbol("b", m.GlobalScope, 1, nothing),
                    m.MonkeySymbol("c", m.FreeScope, 0, nothing),
                    m.MonkeySymbol("d", m.FreeScope, 1, nothing),
                    m.MonkeySymbol("e", m.LocalScope, 0, nothing),
                    m.MonkeySymbol("f", m.LocalScope, 1, nothing),
                ],
                [
                    m.MonkeySymbol("c", m.LocalScope, 0, nothing),
                    m.MonkeySymbol("d", m.LocalScope, 1, nothing),
                ],
            ),
        ]
            for sym in expected_symbols
                @test m.resolve(table, sym.name)[1] == sym
            end

            @test length(table.free_symbols) == length(expected_free_symbols)

            for (actual, expected) in zip(table.free_symbols, expected_free_symbols)
                @test actual == expected
            end
        end
    end

    @testset "Unresolvable Free Variables" begin
        g = m.SymbolTable()
        m.define!(g, "a")

        l1 = m.SymbolTable(; outer = g)
        m.define!(l1, "c")

        l2 = m.SymbolTable(; outer = l1)
        m.define!(l2, "e")
        m.define!(l2, "f")

        expected = [
            m.MonkeySymbol("a", m.GlobalScope, 0, nothing),
            m.MonkeySymbol("c", m.FreeScope, 0, nothing),
            m.MonkeySymbol("e", m.LocalScope, 0, nothing),
            m.MonkeySymbol("f", m.LocalScope, 1, nothing),
        ]

        for sym in expected
            @test m.resolve(l2, sym.name)[1] == sym
        end

        for name in ["b", "d"]
            @test isnothing(m.resolve(l2, name)[1])
        end
    end

    @testset "Define And Resolve Function Name" begin
        g = m.SymbolTable()
        m.define_function!(g, "a")

        @test m.resolve(g, "a")[1] == m.MonkeySymbol("a", m.FunctionScope, 0, nothing)
    end

    @testset "Shadowing Function Name" begin
        g = m.SymbolTable()
        m.define_function!(g, "a")
        m.define!(g, "a")

        @test m.resolve(g, "a")[1] == m.MonkeySymbol("a", m.GlobalScope, 0, nothing)
    end
end


text:
module TestAcquisition
using Test
using LinearAlgebra
using GaussianDistributions
using CovarianceFunctions
const Kernel = CovarianceFunctions
using SARA: ucb, inner_sampling, random_sampling, uncertainty_sampling,
            integrated_uncertainty_sampling

@testset "acquisition" begin
    l = 1/2
    k = Kernel.Lengthscale(Kernel.EQ(), l)
    G = Gaussian(sin, k)
    n = 256

    @testset "random sampling" begin
        x = randn(n)
        @test random_sampling(G, x) isa Real
        @test random_sampling(G, x) in x
    end

    @testset "upper confidence bound" begin
        # test on unconditioned GP
        x = range(0, π, length = n)
        α = 3.
        xi = ucb(G, x, α)
        @test isapprox(xi, π/2, atol = 5e-2)

        # create 2d synthetic data and condition
        l = 1/2
        k = Kernel.Lengthscale(Kernel.EQ(), l)
        G = Gaussian(k)
        f(x) = sum(sin, x)
        x = [randn(2) for _ in 1:n]
        σ = .01
        y = @. f(x) + σ*randn()
        C = G | (x, y, σ^2)

        # test ucb sampling on conditioned process
        ns = 1024
        xs = [randn(2) for _ in 1:ns]

        x0 = ucb(C, xs, 0.) # without uncertainty term, chooses point close to optimum
        @test isapprox(f(x0), 2, atol = 3e-1)

        xi = ucb(C, xs, α) # with uncertainty term,
        @test f(x0) > f(xi) # trades off value with uncertainty
        @test f(xi) > 1 # still chooses xi with moderately large value
    end

    @testset "integrated uncertainty sampling" begin
        l = .1
        k = Kernel.Lengthscale(Kernel.EQ(), l)
        G = Gaussian(k)

        # test on unconditioned GP
        n = 32
        x = range(0, π, length = n)
        σ = 1e-2
        y = @. sin(x) + σ * randn()
        C = G | (x, y, σ^2)

        m = 256
        xs = range(0, π, length = m)

        xi = integrated_uncertainty_sampling(G, xs, σ)
        xk = uncertainty_sampling(G, xs)

        xi = integrated_uncertainty_sampling(C, xs, σ)
        xk = uncertainty_sampling(C, xs)
    end

    @testset "SARA inner loop sampling" begin
        l = .1
        k = Kernel.Lengthscale(Kernel.EQ(), l)
        G = Gaussian(sin, k)

        # test on unconditioned GP
        n = 32
        x = range(0, π, length = n)
        σ = 1e-2
        y = @. sin(x) + σ * randn()
        C = G | (x, y, σ^2)
        m = 256
        xs = range(0, π, length = m)

        xi = inner_sampling(G, xs, σ)
        xi = inner_sampling(C, xs, σ)
        # @test isapprox(xi, π/2, atol = 5e-2)
    end
end

end


text:
# little-endian
@inline bincat(x::UInt8, y::UInt8) = convert(UInt16, x) + (convert(UInt16, y) << 8)

hexstring(n::Unsigned) = sprint(show, n)

# as simple as this seems in retrospect, it was confusing to figure out, hence the shortcut
@inline overflow(a::T, b::T) where T<:Unsigned = b > typemax(T) - a


text:
# In this tutorial, we will look at how to
# - Evaluate `CellFields` at arbitrary points
# - Interpolate finite element functions defined on different
# triangulations. We will consider examples for
#    - Lagrangian finite element spaces
#    - Raviart Thomas finite element spaces
#    - Vector-Valued Spaces
#    - Multifield finite element spaces

# ## Problem Statement
# Let $\mathcal{T}_1$ and $\mathcal{T}_2$ be two triangulations of a
# domain $\Omega$. Let $V_i$ be the finite element space defined on
# the triangulation $\mathcal{T}_i$ for $i=1,2$. Let $f_h \in V_1$. The
# interpolation problem is to find $g_h \in V_2$ such that
#
# ```math
# dof_k^{V_2}(g_h) = dof_k^{V_2}(f_h),\quad \forall k \in
# \{1,\dots,N_{dof}^{V_2}\}
# ```

# ## Setup
# For the purpose of this tutorial we require `Test`, `Gridap` along with the
# following submodules of `Gridap`

using Test
using Gridap
using Gridap.CellData
using Gridap.Visualization

# We now create a computational domain on the unit square $[0,1]^2$ consisting
# of 5 cells per direction

domain = (0,1,0,1)
partition = (5,5)
𝒯₁ = CartesianDiscreteModel(domain, partition)

# ## Background
# `Gridap` offers the feature to evaluate functions at arbitrary
# points in the domain. This will be shown in the next
# section. Interpolation then takes advantage of this feature to
# obtain the `FEFunction` in the new space from the old one by
# evaluating the appropriate degrees of freedom. Interpolation works
# using the composite type `Interpolable` to tell `Gridap` that the
# argument can be interpolated between triangulations.

# ## Interpolating between Lagrangian FE Spaces

# Let us define the infinite dimensional function

f(x) = x[1] + x[2]

# This function will be interpolated to the source `FESpace`
# $V_1$. The space can be built using

reffe₁ = ReferenceFE(lagrangian, Float64, 1)
V₁ = FESpace(𝒯₁, reffe₁)

# Finally to build the function $f_h$, we do

fₕ = interpolate_everywhere(f,V₁)

# To construct arbitrary points in the domain, we use `Random` package:

using Random
pt = Point(rand(2))
pts = [Point(rand(2)) for i in 1:3]

# The finite element function $f_h$ can be evaluated at arbitrary points (or
# array of points) by

fₕ(pt), fₕ.(pts)

# We can also check our results using

@test fₕ(pt) ≈ f(pt)
@test fₕ.(pts) ≈ f.(pts)

# Now let us define the new triangulation $\mathcal{T}_2$ of
# $\Omega$. We build the new triangulation using a partition of 20 cells per
# direction. The map can be passed as an argument to
# `CartesianDiscreteModel` to define the position of the vertices in
# the new mesh.

partition = (20,20)
𝒯₂ = CartesianDiscreteModel(domain,partition)

# As before, we define the new `FESpace` consisting of second order
# elements

reffe₂ = ReferenceFE(lagrangian, Float64, 2)
V₂ = FESpace(𝒯₂, reffe₂)

# Now we interpolate $f_h$ onto $V_2$ to obtain the new function
# $g_h$. The first step is to create the `Interpolable` version of
# $f_h$.

ifₕ = Interpolable(fₕ)

# Then to obtain $g_h$, we dispatch `ifₕ` and the new `FESpace` $V_2$
# to the `interpolate_everywhere` method of `Gridap`.

gₕ = interpolate_everywhere(ifₕ, V₂)

# We can also use
# `interpolate` if interpolating only on the free dofs or
# `interpolate_dirichlet` if interpolating the Dirichlet dofs of the
# `FESpace`.

ḡₕ = interpolate(ifₕ, V₂)

# The finite element function $\bar{g}_h$ is the same as $g_h$ in this
# example since all the dofs are free.

@test gₕ.cell_dof_values ==  ḡₕ.cell_dof_values

# Now we obtain a finite element function using `interpolate_dirichlet`

g̃ₕ = interpolate_dirichlet(ifₕ, V₂)

# Now $\tilde{g}_h$ will be equal to 0 since there are
# no Dirichlet nodes defined in the `FESpace`. We can check by running

g̃ₕ.cell_dof_values

# Like earlier we can check our results for `gₕ`:

@test fₕ(pt) ≈ gₕ(pt) ≈ f(pt)
@test fₕ.(pts) ≈ gₕ.(pts) ≈ f.(pts)

# We can visualize the results using Paraview

writevtk(get_triangulation(fₕ), "source", cellfields=["fₕ"=>fₕ])
writevtk(get_triangulation(gₕ), "target", cellfields=["gₕ"=>gₕ])

# which produces the following output

# ![Target](../assets/interpolation_fe/source_and_target.png)

# ## Interpolating between Raviart-Thomas FESpaces

# The procedure is identical to Lagrangian finite element spaces, as
# discussed in the previous section. The extra thing here is that
# functions in Raviart-Thomas spaces are vector-valued. The degrees of
# freedom of the RT spaces are fluxes of the function across the edge
# of the element. Refer to the
# [tutorial](https://gridap.github.io/Tutorials/dev/pages/t007_darcy/)
# on Darcy equation with RT for more information on the RT
# elements.

# Assuming a function

f(x) = VectorValue([x[1], x[2]])

# on the domain, we build the associated finite dimensional version
# $f_h \in V_1$.

reffe₁ = ReferenceFE(raviart_thomas, Float64, 1) # RT space of order 1
V₁ = FESpace(𝒯₁, reffe₁)
fₕ = interpolate_everywhere(f, V₁)

# As before, we can evaluate the RT function on any arbitrary point in
# the domain.

fₕ(pt), fₕ.(pts)

# Constructing the target RT space and building the `Interpolable`
# object,

reffe₂ = ReferenceFE(raviart_thomas, Float64, 1) # RT space of order 1
V₂ = FESpace(𝒯₂, reffe₂)
ifₕ = Interpolable(fₕ)

# we can construct the new `FEFunction` $g_h \in V_2$ from $f_h$

gₕ = interpolate_everywhere(ifₕ, V₂)

# Like earlier we can check our results

@test gₕ(pt) ≈ f(pt) ≈ fₕ(pt)

# ## Interpolating vector-valued functions

# We can also interpolate vector-valued functions across
# triangulations. First, we define a vector-valued function on a
# two-dimensional mesh.

f(x) = VectorValue([x[1], x[1]+x[2]])

# We then create a vector-valued reference element containing linear
# elements along with the source finite element space $V_1$.

reffe₁ = ReferenceFE(lagrangian, VectorValue{2,Float64}, 1)
V₁ = FESpace(𝒯₁, reffe₁)
fₕ = interpolate_everywhere(f, V₁)

# The target finite element space $V_2$ can be defined in a similar manner.

reffe₂ = ReferenceFE(lagrangian, VectorValue{2,Float64}, 2)
V₂ = FESpace(𝒯₂, reffe₂)

# The rest of the process is similar to the previous sections, i.e.,
# define the `Interpolable` version of $f_h$ and use
# `interpolate_everywhere` to find $g_h \in V₂$.

ifₕ = Interpolable(fₕ)
gₕ = interpolate_everywhere(ifₕ, V₂)

# We can then check the results

@test gₕ(pt) ≈ f(pt) ≈ fₕ(pt)

# ## Interpolating Multi-field Functions

# Similarly, it is possible to interpolate between multi-field finite element
# functions. First, we define the components $h_1(x), h_2(x)$ of a
# multi-field function $h(x)$ as follows.

h₁(x) = x[1]+x[2]
h₂(x) = x[1]

# Next we create a Lagrangian finite element space containing linear
# elements.

reffe₁ = ReferenceFE(lagrangian, Float64, 1)
V₁ = FESpace(𝒯₁, reffe₁)

# Next we create a `MultiFieldFESpace` $V_1 \times V_1$ and
# interpolate the function $h(x)$ to the source space $V_1$.

V₁xV₁ = MultiFieldFESpace([V₁,V₁])
fₕ = interpolate_everywhere([h₁, h₂], V₁xV₁)

# Similarly, the target multi-field finite element space is created
# using $\Omega_2$.

reffe₂ = ReferenceFE(lagrangian, Float64, 2)
V₂ = FESpace(𝒯₂, reffe₂)
V₂xV₂ = MultiFieldFESpace([V₂,V₂])

# Now, to find $g_h \in V_2 \times V_2$, we first extract the components of
# $f_h$ and obtain the `Interpolable` version of the components.

fₕ¹, fₕ² = fₕ
ifₕ¹ = Interpolable(fₕ¹)
ifₕ² = Interpolable(fₕ²)

# We can then use `interpolate_everywhere` on the `Interpolable`
# version of the components and obtain $g_h \in V_2 \times V_2$ as
# follows.

gₕ = interpolate_everywhere([ifₕ¹,ifₕ²], V₂xV₂)

# We can then check the results of the interpolation, component-wise.

gₕ¹, gₕ² = gₕ
@test fₕ¹(pt) ≈ gₕ¹(pt)
@test fₕ²(pt) ≈ gₕ²(pt)

# ## Acknowledgements

# Gridap contributors acknowledge support received from Google,
# Inc. through the Google Summer of Code 2021 project [A fast finite
# element interpolator in
# Gridap.jl](https://summerofcode.withgoogle.com/projects/#6175012823760896).


text:
"""
Project `x` onto sparsity set with `k` non-zero elements.
Assumes `idx` enters as a vector of indices into `x`.
"""
function project_l0_ball!(x, idx, k, buffer)
    n = length(x)
    # do nothing if k > length(x)
    if k ≥ n return x end

    # fill with zeros if k ≤ 0
    if k ≤ 0 return fill!(x, 0) end

    # otherwise, find the spliting element
    search_by_top_k = k < n-k+1
    if search_by_top_k
        _k = k
        pivot = l0_search_partialsort!(idx, x, _k, true)
    else
        _k = n-k+1
        pivot = l0_search_partialsort!(idx, x, _k, false)
    end

    # preserve the top k elements
    p = abs(pivot)
    nonzero_count = 0
    @inbounds for i in eachindex(x)
        if x[i] == 0 continue end
        if abs(x[i]) < p
            x[i] = 0
        else
            nonzero_count += 1
        end
    end

    # resolve ties
    if nonzero_count > k
        number_to_drop = nonzero_count - k
        _buffer_ = view(buffer, 1:number_to_drop)
        _indexes_ = findall(!iszero, x)
        sample!(_indexes_, _buffer_, replace=false)
        @inbounds for i in _buffer_
            x[i] = 0
        end
    end

    return x
end

function project_l0_ball!(X::AbstractMatrix, idx, scores, k, buffer; by::Union{Val{:row}, Val{:col}}=Val(:row))
    # determine structure of sparsity
    if by isa Val{:row}
        n = size(X, 1)
        itr = axes(X, 1)
        itr2 = eachrow(X)
        f = i -> norm(view(X, i, :))
    elseif by isa Val{:col}
        n = size(X, 2)
        itr = axes(X, 2)
        itr2 = eachcol(X)
        f = i -> norm(view(X, :, i))
    else
        error("uncrecognized option `by=$(by)`.")
    end

    # do nothing if k > length(x)
    if k ≥ n return X end

    # fill with zeros if k ≤ 0
    if k ≤ 0 return fill!(X, 0) end

    # otherwise, map rows to a score used in ranking and find the spliting element
    map!(f, scores, itr)
    search_by_top_k = k < n-k+1
    if search_by_top_k
        _k = k
        pivot = l0_search_partialsort!(idx, scores, _k, true)
    else
        _k = n-k+1
        pivot = l0_search_partialsort!(idx, scores, _k, false)
    end

    # preserve the top k elements
    p = abs(pivot)
    nonzero_count = 0
    @inbounds for (i, xᵢ) in enumerate(itr2)
        if scores[i] == 0 continue end

        # row is not in the top k
        if scores[i] < p
            fill!(xᵢ, 0)
            scores[i] = 0
        else # row is in the top k
            nonzero_count += 1
        end
    end

    # resolve ties
    if nonzero_count > k
        number_to_drop = nonzero_count - k
        _buffer_ = view(buffer, 1:number_to_drop)
        _indexes_ = findall(!iszero, scores)
        sample!(_indexes_, _buffer_, replace=false)
        @inbounds for i in _buffer_
            if by isa Val{:row}
                fill!(view(X, i, :), 0)
            elseif by isa Val{:col}
                fill!(view(X, :, i), 0)
            end
        end
    end

    return X
end

"""
Search `x` for the pivot that splits the vector into the `k`-largest elements in magnitude.

The search preserves signs and returns `x[k]` after partially sorting `x`.
"""
function l0_search_partialsort!(idx, x, k, rev::Bool)
    #
    # Based on https://github.com/JuliaLang/julia/blob/788b2c77c10c2160f4794a4d4b6b81a95a90940c/base/sort.jl#L863
    # This eliminates a mysterious allocation of ~48 bytes per call for
    #   sortperm!(idx, x, alg=algorithm, lt=isless, by=abs, rev=true, initialized=false)
    # where algorithm = PartialQuickSort(lo:hi)
    # Savings are small in terms of performance but add up for CV code.
    #
    lo = k
    hi = k+1

    # Order arguments
    lt  = isless
    by  = abs
    # rev = true
    o = Base.Order.Forward
    order = Base.Order.Perm(Base.Sort.ord(lt, by, rev, o), x)

    # sort!(idx, lo, hi, PartialQuickSort(k), order)
    Base.Sort.Float.fpsort!(idx, PartialQuickSort(lo:hi), order)

    return x[idx[k]]
end

struct L0Projection <: Function
    idx::Vector{Int}
    buffer::Vector{Int}

    function L0Projection(n::Int)
        new(collect(1:n), Vector{Int}(undef, n))
    end
end

function (P::L0Projection)(x, k)
    project_l0_ball!(x, P.idx, k, P.buffer)
end

struct StructuredL0Projection <: Function
    idx::Vector{Int}
    buffer::Vector{Int}
    scores::Vector{Float64}

    function StructuredL0Projection(n::Int)
        new(collect(1:n), Vector{Int}(undef, n), Vector{Float64}(undef, n))
    end
end

function (P::StructuredL0Projection)(X::AbstractMatrix, k)
    project_l0_ball!(X, P.idx, P.scores, k, P.buffer, by=Val(:row))
end


text:
function merge_data(idx1::Integer, idx2::Integer, a::Int64, b::Int64)
    return a*b
end
function merge_edges(idx1::Integer, idx2::Integer, d::Float64)
    return 1.0
end
    r = Graph{Integer, Float64}()
    for i = 1:8
        add_vertex!(r, i )
    end
    connect!(r, 1 ,2, 1.0)
    connect!(r, 1 ,3, 1.0)
    connect!(r, 2 ,3, 1.0)
    connect!(r, 1 ,4, 1.0)
    connect!(r, 1 ,5, 1.0)
    connect!(r, 1 ,8, 1.0)
    connect!(r, 4 ,6, 1.0)
    connect!(r, 5 ,7, 1.0)

    merge_vertices!(r, 4 ,3)

    [2;5;8;9] == sort([edge.node.index for edge in r.nodes[1].neighbors])
    [1;9] == sort([edge.node.index for edge in r.nodes[2].neighbors])
    [1;2;6]== sort([edge.node.index for edge in r.nodes[9].neighbors])
    [1;7]== sort([edge.node.index for edge in r.nodes[5].neighbors])
    [1]== sort([edge.node.index for edge in r.nodes[8].neighbors])
    [5]== sort([edge.node.index for edge in r.nodes[7].neighbors])
    [9]== sort([edge.node.index for edge in r.nodes[6].neighbors])


text:
mutable struct MosekLinearQuadraticModel <: MathProgBase.AbstractLinearQuadraticModel
    task :: Mosek.MSKtask

    binvarflags:: Vector{Bool}
    # NOTE: bkx/blx/bux are the bound on the variable in the
    # continuous problem. Setting a variable to :Bin, will not change
    # these. Setting :Cont or :Int on a :Bin variable will effectively
    # revert to the original bounds.
    numvar     :: Int
    numcon     :: Int
    bkx        :: Vector{Mosek.Boundkey}
    blx        :: Vector{Float64}
    bux        :: Vector{Float64}

    bkc        :: Vector{Mosek.Boundkey}
    blc        :: Vector{Float64}
    buc        :: Vector{Float64}
    lincon     :: Vector{Int32}
    quadcon    :: Vector{Int32}

    lasttrm    :: Mosek.Rescode

    options
end

mutable struct MosekNonlinearModel <: MathProgBase.AbstractLinearQuadraticModel
    m :: MosekLinearQuadraticModel
end

function MathProgBase.LinearQuadraticModel(s::Mosek.MosekSolver)
  r = MosekLinearQuadraticModel(Mosek.maketask(),
                                Bool[],
                                0,
                                0,
                                Int32[],
                                Float64[],
                                Float64[],

                                Int32[],
                                Float64[],
                                Float64[],

                                Int32[],
                                Int32[],

                                Mosek.MSK_RES_OK,

                                s.options)
    loadoptions!(r)
    r
end
##############################################################
## Linear
##############################################################

function MathProgBase.loadproblem!(m::MosekLinearQuadraticModel,
                                   A,
                                   collb :: Vector{T1},
                                   colub :: Vector{T2},
                                   obj   :: Vector{T3},
                                   rowlb :: Vector{T4},
                                   rowub :: Vector{T5},
                                   sense:: Symbol) where {T1,T2,T3,T4,T5}
    MathProgBase.loadproblem!(m,
                              convert(SparseMatrixCSC{Float64,Int},A),
                              convert(Vector{Float64},collb),
                              convert(Vector{Float64},colub),
                              convert(Vector{Float64},obj),
                              convert(Vector{Float64},rowlb),
                              convert(Vector{Float64},rowub),
                              sense)
end

function MathProgBase.loadproblem!(m::MosekLinearQuadraticModel,
                                   A::     SparseMatrixCSC{Float64,Int},
                                   collb:: Vector{Float64},
                                   colub:: Vector{Float64},
                                   obj::   Vector{Float64},
                                   rowlb:: Vector{Float64},
                                   rowub:: Vector{Float64},
                                   sense:: Symbol)
    Mosek.resizetask(m.task,0,0,0,0,0);

    nrows,ncols = size(A)
    if ncols != length(collb) ||
        ncols != length(colub) ||
        ncols != size(obj,1)   ||
        nrows != length(rowlb) ||
        nrows != length(rowub) ||
        ncols != length(obj)

        throw(MosekMathProgModelError("Inconsistent data dimensions"))
    end

    Mosek.appendvars(m.task,ncols)
    Mosek.appendcons(m.task,nrows)

    (m.bkx,m.blx,m.bux) = makebounds(collb,colub)

    (m.bkc,m.blc,m.buc) = makebounds(rowlb,rowub)

    m.numvar = length(m.bkx)
    m.numcon = length(m.bkc)
    m.lincon = Int32[1:nrows;]
    m.quadcon = Int32[]
    m.binvarflags = fill(false,m.numvar)

    # input coefficients
    Mosek.putclist(m.task, Int32[1:ncols;], obj)
    Mosek.putacolslice(m.task, Int32(1), Int32(ncols+1), A.colptr[1:ncols], A.colptr[2:ncols+1], A.rowval, A.nzval)
    MathProgBase.setsense!(m, sense)

    # input bounds
    Mosek.putvarboundslice(m.task, Int32(1), Int32(ncols+1), m.bkx, m.blx, m.bux)
    Mosek.putconboundslice(m.task, Int32(1), Int32(nrows+1), m.bkc, m.blc, m.buc)
    m
end

function MathProgBase.loadproblem!(m::MosekLinearQuadraticModel,
                                   filename:: AbstractString)
    tmptask = Mosek.maketask()
    try
        readdata(tmptask, filename)

        if  Mosek.getnumcone(tmptask) > 0 ||
            Mosek.getnumbarvar(tmptask) > 0
            throw(MosekMathProgModelError("Not a linear/quadratic model"))
        end

        numcon = Mosek.getnumcon(tmptask)
        numvar = Mosek.getnumvar(tmptask)

        bkx,blx,bux = Mosek.getvarboundslice(tmptask,Int32(1),Int32(numvar+1))
        bkc,blc,buc = Mosek.getconboundslice(tmptask,Int32(1),Int32(numcon+1))

        vts = Mosek.getvartypelist(tmptask,Int32[1:numvar])
        binflags = Bool[ (vts[i] == Mosek.MSK_VARIABLE_TYPE_INT &&
                          bkx[i] == Mosek.MSK_BK_RA &&
                          abs(blx[i]) < 1e-8 &&
                          abs(bux[i]-1.0) < 1e-8)
                        for i in 1:numvar ]

        lincon = findall(i -> Mosek.getnumqconknz(m.task,i) == 0, 1:numcon)
        quadcon = findall(i -> Mosek.getnumqconknz(m.task,i) > 0, 1:numcon)

        Mosek.deletetask(m.task)

        m.task = tmptask
        m.binvaflags = binflags
        m.bkx = bkx
        m.blx = blx
        m.bux = bux
        m.bkc = bkc
        m.blc = blc
        m.buc = buc
        m.numvar = lenght(bkx)
        m.numcon = length(bkc)
        m.lincon = lincon
        m.quadcon = quadcon
    catch
        Mosek.deletetask(tmptask)
        rethrow()
    end

    m
end

function loadoptions!(m::MosekLinearQuadraticModel)
  loadoptions_internal!(m.task, m.options)
end

function MathProgBase.writeproblem(m::MosekLinearQuadraticModel, filename::AbstractString)
    Mosek.putintparam(m.task,Mosek.MSK_IPAR_OPF_WRITE_SOLUTIONS, Mosek.MSK_ON)
    Mosek.writedata(m.task,filename)
end

function MathProgBase.getvarLB(m::MosekLinearQuadraticModel)
    bk,bl,bu = Mosek.getvarboundslice(m.task,Int32(1),Int32(m.numvar+1))
    for i in 1:length(bk)
        if bk[i] == Mosek.MSK_BK_FR || bk[i] == Mosek.MSK_BK_UP
            bl[i] = -Inf
        end
    end
    bl
end

function MathProgBase.getvarUB(m::MosekLinearQuadraticModel)
    bk,bl,bu = Mosek.getvarboundslice(m.task,Int32(1),Int32(m.numvar+1))
    for i in 1:length(bk)
        if bk[i] == Mosek.MSK_BK_FR || bk[i] == Mosek.MSK_BK_LO
            bu[i] = Inf
        end
    end
    bu
end

function MathProgBase.getconstrLB(m::MosekLinearQuadraticModel)
    bk,bl,bu = Mosek.getconboundslice(m.task,1,m.numcon+1)
    for i in 1:length(bk)
        if bk[i] == Mosek.MSK_BK_FR || bk[i] == Mosek.MSK_BK_UP
            bl[i] = -Inf
        end
    end
    bl
end

function MathProgBase.getconstrUB(m::MosekLinearQuadraticModel)
    bk,bl,bu = Mosek.getconboundslice(m.task,1,m.numcon+1)
    for i in 1:length(bk)
        if bk[i] == Mosek.MSK_BK_FR || bk[i] == Mosek.MSK_BK_LO
            bu[i] = Inf
        end
    end
    bu
end

MathProgBase.setvarLB!(m::MosekLinearQuadraticModel, bnd::Vector{T}) where {T} = MathProgBase.setvarLB!(m,convert(Vector{Float64},bnd))
function MathProgBase.setvarLB!(m::MosekLinearQuadraticModel, bnd::Vector{Float64})
    n = min(length(bnd),m.numvar)

    for i in 1:n
        if bnd[i] > -Inf
            if m.bux[i] < Inf
                if abs(bnd[i]-m.bux[i]) < 1e-8
                    m.bkx[i] = Mosek.MSK_BK_FX
                    m.blx[i] = m.bux[i]
                    m.bux[i] = m.bux[i]
                else
                    m.bkx[i] = Mosek.MSK_BK_RA
                    m.blx[i] = bnd[i]
                end
            else # bux[i] == Inf
                m.bkx[i] = Mosek.MSK_BK_LO
                m.blx[i] = bnd[i]
            end
        else # bnd[i] == -Inf
            if m.bux[i] < Inf
                m.bkx[i] = Mosek.MSK_BK_UP
            else
                m.bkx[i] = Mosek.MSK_BK_FR
            end
            m.blx[i] = -Inf
        end
    end

    Mosek.putvarboundslice(m.task,1,n+1,m.bkx,m.blx,m.bux)
    if any(m.binvarflags)
        idxs = convert(Vector{Int32},findall(v->v, m.binvarflags))
        bkx = Mosek.Boundkey[ Mosek.MSK_BK_RA for i in 1:length(idxs)]
        blx = Float64[ max(m.blx[i],0.0) for i in idxs ]
        bux = Float64[ min(m.bux[i],1.0) for i in idxs ]

        Mosek.putvarboundlist(m.task,idxs, bkx,blx,bux)
    end

    nothing
end

MathProgBase.setvarUB!(m::MosekLinearQuadraticModel, bnd::Vector{T}) where {T} = MathProgBase.setvarUB!(m,convert(Vector{Float64},bnd))
function MathProgBase.setvarUB!(m::MosekLinearQuadraticModel, bnd::Vector{Float64})
    n = min(length(bnd),m.numvar)

    for i in 1:n
        if bnd[i] < Inf
            if m.blx[i] > -Inf
                if abs(bnd[i]-m.blx[i]) < 1e-8
                    m.bkx[i] = Mosek.MSK_BK_FX
                    m.blx[i] = m.blx[i]
                    m.bux[i] = m.blx[i]
                else
                    m.bkx[i] = Mosek.MSK_BK_RA
                    m.bux[i] = bnd[i]
                end
            else # blx[i] == -Inf
                m.bkx[i] = Mosek.MSK_BK_UP
                m.bux[i] = bnd[i]
            end
        else # bnd[i] == Inf
            if m.blx[i] > -Inf
                m.bkx[i] = Mosek.MSK_BK_LO
            else
                m.bkx[i] = Mosek.MSK_BK_FR
            end
            m.bux[i] = Inf
        end
    end

    Mosek.putvarboundslice(m.task,1,n+1,m.bkx,m.blx,m.bux)
    if any(m.binvarflags)
        idxs = convert(Vector{Int32},findall(v->v, m.binvarflags))
        bkx = Mosek.Boundkey[ Mosek.MSK_BK_RA for i in 1:length(idxs)]
        blx = Float64[ max(m.blx[i],0.0) for i in idxs ]
        bux = Float64[ min(m.bux[i],1.0) for i in idxs ]

        Mosek.putvarboundlist(m.task,idxs, bkx,blx,bux)
    end

    nothing
end

MathProgBase.setconstrLB!(m::MosekLinearQuadraticModel, bnd::Vector{T}) where {T} = MathProgBase.setconstrLB!(m,convert(Vector{Float64},bnd))
function MathProgBase.setconstrLB!(m::MosekLinearQuadraticModel, bnd::Vector{Float64})
    n = min(length(bnd),length(m.lincon))

    for i in 1:n
        bk = m.bkc[m.lincon[i]]
        if bnd[i] > -Inf
            if bk == Mosek.MSK_BK_UP || bk == Mosek.MSK_BK_RA || bk == Mosek.MSK_BK_FX
                if abs(bnd[i]-m.buc[m.lincon[i]]) < 1e-8
                    m.bkc[m.lincon[i]] = Mosek.MSK_BK_FX
                    m.blc[m.lincon[i]] = m.buc[m.lincon[i]]
                    m.buc[m.lincon[i]] = m.buc[m.lincon[i]]
                else
                    m.bkc[m.lincon[i]] = Mosek.MSK_BK_RA
                    m.blc[m.lincon[i]] = bnd[i]
                end
            else # buc[m.lincon[i]] == Inf
                m.bkc[m.lincon[i]] = Mosek.MSK_BK_LO
                m.blc[m.lincon[i]] = bnd[i]
            end
        else # bnd[m.lincon[i]] == -Inf
            if bk == Mosek.MSK_BK_UP || bk == Mosek.MSK_BK_RA || bk == Mosek.MSK_BK_FX
                m.bkc[m.lincon[i]] = Mosek.MSK_BK_UP
            else
                m.bkc[m.lincon[i]] = Mosek.MSK_BK_FR
            end
            m.blc[m.lincon[i]] = -Inf
        end
    end

    Mosek.putconboundslice(m.task,Int32(1),Int32(m.numcon+1),m.bkc,m.blc,m.buc)

    nothing
end

MathProgBase.setconstrUB!(m::MosekLinearQuadraticModel, bnd::Vector{T}) where {T} = MathProgBase.setconstrUB!(m,convert(Vector{Float64},bnd))
function MathProgBase.setconstrUB!(m::MosekLinearQuadraticModel, bnd::Vector{Float64})
    n = min(length(bnd),m.numcon)

    for i in 1:n
        bk = m.bkc[m.lincon[i]]
        if bnd[i] < Inf
            if bk == Mosek.MSK_BK_LO || bk == Mosek.MSK_BK_RA || bk == Mosek.MSK_BK_FX
                if abs(bnd[i]-m.blc[m.lincon[i]]) < 1e-8
                    m.bkc[m.lincon[i]] = Mosek.MSK_BK_FX
                    m.buc[m.lincon[i]] = m.blc[m.lincon[i]]
                else
                    m.bkc[m.lincon[i]] = Mosek.MSK_BK_RA
                    m.buc[m.lincon[i]] = bnd[i]
                end
            else # blc[i] == -Inf
                m.bkc[m.lincon[i]] = Mosek.MSK_BK_UP
                m.buc[m.lincon[i]] = bnd[i]
            end
        else # bnd[i] == Inf
            if bk == Mosek.MSK_BK_LO || bk == Mosek.MSK_BK_RA || bk == Mosek.MSK_BK_FX
                m.bkc[m.lincon[i]] = Mosek.MSK_BK_LO
            else
                m.bkc[m.lincon[i]] = Mosek.MSK_BK_FR
            end
            m.buc[m.lincon[i]] = Inf
        end
    end
    Mosek.putconboundslice(m.task,Int32(1),Int32(m.numcon+1),m.bkc,m.blc,m.buc)

    nothing
end

function MathProgBase.getobj(m::MosekLinearQuadraticModel)
    Mosek.getcslice(m.task,1,m.numvar+1)
end

MathProgBase.setobj!(m::MosekLinearQuadraticModel, c::Vector{T}) where {T} = MathProgBase.setobj!(m,convert(Vector{Float64},c))
function MathProgBase.setobj!(m::MosekLinearQuadraticModel, c :: Vector{Float64})
    n = min(length(c),m.numvar)
    Mosek.putclist(m.task,Int32[1:n;],c[1:n])
end

function MathProgBase.getconstrmatrix(m::MosekLinearQuadraticModel)
    (ptrb,ptre,asubj,aval) = Mosek.getaslice(m.task, Mosek.MSK_ACC_CON,1,m.numcon+1)

    numnz = length(asubj)
    asubi = Vector{Int32}(undef,numnz)

    for i in 1:m.numcon
        asubi[ptrb[i]:ptre[i]-1] .= Int32(i)
    end

    sparse(asubi,asubj,aval,m.numcon,m.numvar)
end

MathProgBase.addvar!(m::MosekLinearQuadraticModel, bl, bu, c) = MathProgBase.addvar!(m,convert(Float64,bl),convert(Float64,bu),convert(Float64,c))

function MathProgBase.addvar!(m::MosekLinearQuadraticModel,
                 bl  ::Float64,
                 bu  ::Float64,
                 c   ::Float64)
    m.numvar += 1
    if bl > -Inf
        if bu < Inf
            if abs(bl-bu) < 1e-8
                push!(m.bkx,Mosek.MSK_BK_FX)
                push!(m.blx,bl)
                push!(m.bux,bl)
            else
                push!(m.bkx,Mosek.MSK_BK_RA)
                push!(m.blx,bl)
                push!(m.bux,bu)
            end
        else
            push!(m.bkx,Mosek.MSK_BK_LO)
            push!(m.blx,bl)
            push!(m.bux,Inf)
        end
    else
        if bu < Inf
            push!(m.bkx,Mosek.MSK_BK_UP)
            push!(m.blx,-Inf)
            push!(m.bux,bu)
        else
            push!(m.bkx,Mosek.MSK_BK_FR)
            push!(m.blx,-Inf)
            push!(m.bux,Inf)
        end
    end

    Mosek.appendvars(m.task,1);
    Mosek.putvarbound(m.task,m.numvar,m.bkx[m.numvar],m.blx[m.numvar],m.bux[m.numvar])
    Mosek.putcj(m.task,m.numvar,c)
    push!(m.binvarflags,false)
end

MathProgBase.addvar!(m::MosekLinearQuadraticModel,
                     subi::Vector{T1},
                     val ::Vector{T2},
                     bl  ::T3,
                     bu  ::T4,
                     c   ::T5) where {T1,T2,T3,T4,T5} = MathProgBase.addvar!(m,convert(Vector{Int32},subi),convert(Vector{Float64},val),convert(Float64,bl),convert(Float64,bu),convert(Float64,c))
function MathProgBase.addvar!(m::MosekLinearQuadraticModel,
                              subi::Vector{Int32},
                              val ::Vector{Float64},
                              bl  ::Float64,
                              bu  ::Float64,
                              c   ::Float64)
    MathProgBase.addvar!(m,bl,bu,c)
    Mosek.putacol(m.task,m.numvar,subi,val)
end

MathProgBase.addconstr!(m::MosekLinearQuadraticModel,
                        subj::Vector{T1},
                        val ::Vector{T2},
                        bl  ::T3,
                        bu  ::T4) where {T1,T2,T3,T4} = MathProgBase.addconstr!(m,convert(Vector{Int32},subj),convert(Vector{Float64},val),convert(Float64,bl),convert(Float64,bu))

function MathProgBase.addconstr!(m::MosekLinearQuadraticModel,
                                 subj::Vector{Int32},
                                 val ::Vector{Float64},
                                 bl  ::Float64,
                                 bu  ::Float64)
    m.numcon += 1
    push!(m.lincon,m.numcon)
    if bl > -Inf
        if bu < Inf
            if abs(bl-bu) < 1e-8
                push!(m.bkc,Mosek.MSK_BK_FX)
                push!(m.blc,bl)
                push!(m.buc,bl)
            else
                push!(m.bkc,Mosek.MSK_BK_RA)
                push!(m.blc,bl)
                push!(m.buc,bu)
            end
        else
            push!(m.bkc,Mosek.MSK_BK_LO)
            push!(m.blc,bl)
            push!(m.buc,Inf)
        end
    else
        if bu < Inf
            push!(m.bkc,Mosek.MSK_BK_UP)
            push!(m.blc,-Inf)
            push!(m.buc,bu)
        else
            push!(m.bkc,Mosek.MSK_BK_FR)
            push!(m.blc,-Inf)
            push!(m.buc,Inf)
        end
    end
    Mosek.appendcons(m.task,1);
    Mosek.putconbound(m.task,m.numcon,m.bkc[m.numcon],m.blc[m.numcon],m.buc[m.numcon])
    Mosek.putarow(m.task,m.numcon,subj,val)
end

MathProgBase.numlinconstr(m::MosekLinearQuadraticModel) = length(m.lincon)

MathProgBase.getobjval(m::MosekLinearQuadraticModel) = getobjval(m.task)

function MathProgBase.getsolution(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)

    solsta = Mosek.getsolsta(m.task,sol)

    if solsta in [Mosek.MSK_SOL_STA_OPTIMAL,
                  Mosek.MSK_SOL_STA_PRIM_FEAS,
                  Mosek.MSK_SOL_STA_PRIM_AND_DUAL_FEAS,
                  Mosek.MSK_SOL_STA_INTEGER_OPTIMAL ]
        Mosek.getxx(m.task,sol)
  else
    throw(MosekMathProgModelError("No solution available"))
  end
end

function MathProgBase.getconstrsolution(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)

    Mosek.getxc(m.task,sol)[m.lincon]
end

function MathProgBase.getreducedcosts(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)
    solsta = Mosek.getsolsta(m.task,sol)

    if solsta in [Mosek.MSK_SOL_STA_OPTIMAL,
                  Mosek.MSK_SOL_STA_DUAL_FEAS,
                  Mosek.MSK_SOL_STA_PRIM_AND_DUAL_FEAS ]
        Mosek.getslx(m.task,sol) - Mosek.getsux(m.task,sol)
    else
        throw(MosekMathProgModelError("No solution available"))
    end
end

function MathProgBase.getconstrduals(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)
    if sol == Mosek.MSK_SOL_ITG
        throw(Mosek.MosekMathProgModelError("Solution not available"))
    end
    solsta = Mosek.getsolsta(m.task,sol)

    if solsta in [Mosek.MSK_SOL_STA_OPTIMAL,
                  Mosek.MSK_SOL_STA_DUAL_FEAS,
                  Mosek.MSK_SOL_STA_PRIM_AND_DUAL_FEAS ]
        Mosek.gety(m.task,sol)[m.lincon]
    else
        throw(MosekMathProgModelError("Solution not available"))
    end
end

function MathProgBase.getbasis(m::MosekLinearQuadraticModel)
    if ! Mosek.solutiondef(m.task,Mosek.MSK_SOL_BAS)
        throw(Mosek.MosekMathProgModelError("Basis not available"))
    end
    sol = Mosek.MSK_SOL_BAS

    skx = Mosek.getskx(m.task,sol)
    skc = Mosek.getskc(m.task,sol)

    cbasis = [if     skx[i] == Mosek.MSK_SK_BAS :Basic
              elseif skx[i] == Mosek.MSK_SK_LO  :NonBasicAtLower
              elseif skx[i] == Mosek.MSK_SK_UP  :NonBasicAtUpper
              elseif skx[i] == Mosek.MSK_SK_FX  :NonBasicAtLower # or upper. Doesn't matter.
              else                        :SuperBasic
              end
              for i in 1:m.numvar ]
    rbasis = [if     skc[i] == Mosek.MSK_SK_BAS :Basic
              elseif skc[i] == Mosek.MSK_SK_LO  :NonBasicAtLower
              elseif skc[i] == Mosek.MSK_SK_UP  :NonBasicAtUpper
              elseif skc[i] == Mosek.MSK_SK_FX  :NonBasicAtLower # or upper. Doesn't matter.
              else                        :SuperBasic
              end
              for i in m.lincon]

    cbasis,rbasis
end

function MathProgBase.getunboundedray(m::MosekLinearQuadraticModel)
    soldef = getsoldef(m.task)

    solsta = Mosek.getsolsta(m.task,soldef)
    if solsta in [ Mosek.MSK_SOL_STA_DUAL_INFEAS_CER ]
        Mosek.getxx(m.task,soldef)
    else
        throw(MosekMathProgModelError("No ray available"))
    end
end

function MathProgBase.getinfeasibilityray(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)

    solsta = Mosek.getsolsta(m.task,sol)
    if solsta in [ Mosek.MSK_SOL_STA_PRIM_INFEAS_CER ]
        Mosek.getslc(m.task,sol) - Mosek.getsuc(m.task,sol)
    else
        throw(MosekMathProgModelError("No ray available"))
    end
end

MathProgBase.getrawsolver(m::MosekLinearQuadraticModel) = m.task

MathProgBase.getsimplexiter(m::MosekLinearQuadraticModel) = Mosek.getintinf(m.task,Mosek.MSK_IINF_SIM_PRIMAL_ITER)+Mosek.getintinf(m.task,Mosek.MSK_IINF_SIM_DUAL_ITER)
MathProgBase.getbarrieriter(m::MosekLinearQuadraticModel) = Mosek.getintinf(m.task,Mosek.MSK_IINF_INTPNT_ITER)

MathProgBase.setwarmstart!(m::MosekLinearQuadraticModel, v::Vector{T}) where {T} = MathProgBase.setwarmstart!(m,convert(Vector{Float64},v))
function MathProgBase.setwarmstart!(m::MosekLinearQuadraticModel, v::Vector{Float64})
    n = min(m.numvar,length(v))
    vals = Vector{Float64}(undef,n)
    vals[:] = v[1:n]

    nanidxs = findall(isnan,vals)
    vals[nanidxs] = 0.0

    skx = Mosek.Stakey[ if isnan(vals[i]) Mosek.MSK_SK_UNK else Mosek.MSK_SK_BAS end for i in 1:n ]

    Mosek.putxxslice(m.task,Mosek.MSK_SOL_BAS,1,n+1,vals);
    Mosek.putskxslice(m.task,Mosek.MSK_SOL_BAS,1,n+1,skx);
end

function MathProgBase.optimize!(m::MosekLinearQuadraticModel)
    try
        m.lasttrm = Mosek.optimize(m.task)
        Mosek.solutionsummary(m.task,Mosek.MSK_STREAM_LOG)
    catch err
        if isa(err, Mosek.MosekError)
            m.lasttrm = err.rcode
        else
            rethrow()
        end
    end
end

MathProgBase.status(m::MosekLinearQuadraticModel) = status(m.task,m.lasttrm)

MathProgBase.getobjbound(m::MosekLinearQuadraticModel) = Mosek.getdouinf(m.task,Mosek.MSK_DINF_MIO_OBJ_BOUND)

MathProgBase.getobjgap(m::MosekLinearQuadraticModel) = getobjgap(m.task)

MathProgBase.getsolvetime(m::MosekLinearQuadraticModel) = Mosek.getdouinf(m.task,Mosek.MSK_DINF_OPTIMIZER_TIME)

MathProgBase.getsense(m::MosekLinearQuadraticModel) = getsense(m.task)

MathProgBase.setsense!(m::MosekLinearQuadraticModel,sense) = setsense!(m.task,sense)

function MathProgBase.freemodel!(m::MosekLinearQuadraticModel)
    Mosek.deletetask(m.task)
    nothing
end

MathProgBase.numvar(m::MosekLinearQuadraticModel) = m.numvar
MathProgBase.numconstr(m::MosekLinearQuadraticModel) = length(m.lincon)

function MathProgBase.setvartype!(m::MosekLinearQuadraticModel,vtvec::Vector{Symbol})
    n = min(m.numvar,length(vtvec))
    if n > 0
        vts = Mosek.Variabletype[if     vt == :Cont Mosek.MSK_VAR_TYPE_CONT
                                 elseif vt == :Int  Mosek.MSK_VAR_TYPE_INT
                                 elseif vt == :Bin  Mosek.MSK_VAR_TYPE_INT
                                 else               Mosek.MSK_VAR_TYPE_CONT
                                 end
                    for vt in vtvec[1:n]]
        Mosek.putvartypelist(m.task,Int32[1:n;],vts)
        for i in findall(vt -> vt == :Bin, vtvec[1:n])
            bl = max(m.blx[i],0.0)
            bu = min(m.bux[i],1.0)
            Mosek.putvarbound(m.task,i,Mosek.MSK_BK_RA,bl,bu)
        end

        # for all :Bin vars being changed to :Int or :Cont, restore original bounds
        for i in findall(i -> (vtvec[i] == :Cont || vtvec[i] == :Int) && m.binvarflags[i], 1:n)
            Mosek.putvarbound(m.task,i,m.bkx[i],m.blx[i],m.bux[i])
        end

        for i in 1:n
            m.binvarflags[i] = vtvec[i] == :Bin
        end
    end
end

function MathProgBase.getvartype(m::MosekLinearQuadraticModel)
    mskvt = Mosek.getvartypelist(m.task,Int32[1:m.numvar;])
    [if mskvt[i] == Mosek.MSK_VAR_TYPE_INT
         if m.binvarflags[i]
             :Bin
         else
             :Int
         end
     else
         :Cont
     end
     for i in 1:m.numvar]
end

##############################################################
## Integer Programming
##############################################################

MathProgBase.getnodecount(m::MosekLinearQuadraticModel) = 0

##############################################################
## Quadratic
##############################################################

MathProgBase.numquadconstr(m::MosekLinearQuadraticModel) = length(m.quadcon)

MathProgBase.setquadobj!(m::MosekLinearQuadraticModel,subi,subj,valij) = MathProgBase.setquadobj!(m,convert(Vector{Int32},subi),convert(Vector{Int32},subj),convert(Vector{Float64},valij))

# NOTE on data format: The matrix is specified by inputting only lower
# or upper triangular part. This means that whenever element (i,j) is
# added, (j,i) is implicitly added. Duplicates are added together
function MathProgBase.setquadobj!(m::MosekLinearQuadraticModel,
                                  subi  :: Vector{Int32},
                                  subj  :: Vector{Int32},
                                  valij :: Vector{Float64})
    n = length(subi)
    let qsubi = subi[:],
        qsubj = subj[:]
        for i in 1:n
            if qsubi[i] < qsubj[i]
                tmp = qsubi[i]
                qsubi[i] = qsubj[i]
                qsubj[i] = tmp
            end
        end
        Mosek.putqobj(m.task,qsubi,qsubj,valij)
    end
end

function MathProgBase.addquadconstr!(m :: MosekLinearQuadraticModel,
                                     subj,
                                     valj,
                                     qsubi,
                                     qsubj,
                                     qvalij,
                                     sense  :: Char,
                                     bnd)
    MathProgBase.addquadconstr!(m,
                                convert(Vector{Int32},subj),
                                convert(Vector{Float64},valj),
                                convert(Vector{Int32},qsubi),
                                convert(Vector{Int32},qsubj),
                                convert(Vector{Float64},qvalij),
                                sense,
                                convert(Float64,bnd))
end

function MathProgBase.addquadconstr!(m      :: MosekLinearQuadraticModel,
                                     subj   :: Vector{Int32},
                                     valj   :: Vector{Float64},
                                     qsubi  :: Vector{Int32},
                                     qsubj  :: Vector{Int32},
                                     qvalij :: Vector{Float64},
                                     sense  :: Char,
                                     bnd    :: Float64)

    if     sense == '<'
        push!(m.bkc,Mosek.MSK_BK_UP)
    elseif sense == '>'
        push!(m.bkc,Mosek.MSK_BK_LO)
    else
        throw(MosekMathProgSolverInterface.MosekMathProgModelError("Invalid sense"))
    end

    m.numcon += 1
    push!(m.quadcon,m.numcon)
    push!(m.blc,bnd)
    push!(m.buc,bnd)

    Mosek.appendcons(m.task,1)
    Mosek.putconbound(m.task,m.numcon,m.bkc[m.numcon],m.blc[m.numcon],m.buc[m.numcon])

    let qsubi = qsubi[:],
        qsubj = qsubj[:],
        qval  = qvalij[:]
        for i in 1:length(qsubi)
            if qsubi[i] < qsubj[i]
                t = qsubi[i]
                qsubi[i] = qsubj[i]
                qsubj[i] = t
            elseif  qsubi[i] == qsubj[i]
                qval[i] *= 2
            end
        end
        Mosek.putqconk(m.task,m.numcon,qsubi,qsubj,qval)
    end

    Mosek.putarow(m.task,m.numcon,subj,valj)
end

function MathProgBase.getquadconstrsolution(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)

    Mosek.getxc(sol)[m.quadcon]
end

function MathProgBase.getquadconstrduals(m::MosekLinearQuadraticModel)
    sol = getsoldef(m.task)
    solsta = Mosek.getsolsta(m.task,sol)

    if solsta in [Mosek.MSK_SOL_STA_OPTIMAL,
                  Mosek.MSK_SOL_STA_DUAL_FEAS,
                  Mosek.MSK_SOL_STA_PRIM_AND_DUAL_FEAS ]
        Mosek.gety(m.task)[m.quadcon]
    else
        throw(MosekMathProgModelError("Solution not available"))
    end
end

function MathProgBase.getquadinfeasibilityray(m::MosekLinearQuadraticModel)
    sol = getsoldef(m)
    solsta = getsolsta(m.task,sol)

    s = Mosek.getsux(m.task,sol) - Mosek.getslx(m.task,sol)
    if solsta in [ Mosek.MSK_SOL_STA_PRIM_INFEAS_CER ]
        -s[m.quadcon]
    else
        throw(MosekMathProgModelError("No solution available"))
    end
end

function MathProgBase.getquadconstrRHS(m::MosekLinearQuadraticModel)
    m.blc[m.quadcon]
end

MathProgBase.setquadconstrRHS!(m::MosekLinearQuadraticModel, bnd) = MathProgBase.setquadconstrRHS!(m,convert(Vector{Float64},bnd))
function MathProgBase.setquadconstrRHS!(m::MosekLinearQuadraticModel, bnd::Vector{Float64})
    n = min(length(bnd),length(m.quadcon))
    m.blc[m.quadcon[1:n]] = bnd[1:n]

    Mosek.putconboundlist(m.task,m.quadcon[1:n],m.bkc,m.blc,m.buc)
end

##############################################################
## Nonlinear hasbeen discontinued
#############################################################


text:
using vOptGeneric
using GLPK,GLPKMathProgInterface
m = vModel(solver = GLPKSolverMIP())
# using CPLEX
# m = vModel(solver = CplexSolver())

@variable(m, x[1:6,1:6], Bin)

M = 50
p1 =   [M 4 5 M M M;
        M M 2 1 2 7;
        M M M 5 2 M;
        M M 5 M M 3;
        M M M M M 4;
        M M M M M M]
p2 =   [M 3 1 M M M;
        M M 1 4 2 2;
        M M M 1 7 M;
        M M 1 M M 2;
        M M M M M 2;
        M M M M M M]

@addobjective(m, Min , sum(x[i,j]*p1[i,j] for i=1:6, j=1:6))
@addobjective(m, Min , sum(x[i,j]*p2[i,j] for i=1:6, j=1:6))

@constraint(m, eps, sum(x[1,j] for j = 1:6) == 1)
@constraint(m, sum(x[i,6] for i = 1:6) == 1)
@constraint(m, cstr[i=2:5], sum(x[i,j] for j = 1:6) - sum(x[j,i] for j =1:6) == 0)

solve(m, method=:epsilon, step=0.5)

#get the results and print/plot them
Y_N = getY_N(m)

for n = 1:length(Y_N)
    X = getvalue(x, n)
    for ind in find(X)
        i,j = ind2sub(X,ind)
        print("$i->$j ")
    end
    println("| z = ",Y_N[n])
end

# using PyPlot
# f1 = map(y -> y[1], Y_N)
# f2 = map(y -> y[2], Y_N)
# xlabel("z1")
# ylabel("z2")
# plot(f1,f2,"bx", markersize = "6")
# !isinteractive() && show()


text:
using Glimpse
const Gl = Glimpse
function spherepoint(r, theta, phi)
    return r * Point3f0(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta))
end
nspheres = 20000
radii = Iterators.cycle(range(1.0f0, 200.0f0; length = div(nspheres, 1)))
angs = Iterators.cycle(range(0, 2pi; length = div(nspheres, 1)))
angs2 = Iterators.cycle(range(-pi, pi; length = div(nspheres, 1)))

ps       = [spherepoint(r, a1, a2) for (r, a1, a2, i) in zip(radii, angs, angs2, 1:nspheres)]
cs       = [RGBAf0(0.1, 0.5, 0.9, 0.4), RGBAf0(0.9, 0.1, 0.5, 0.4)]
sph_geom = Gl.PolygonGeometry(Gl.Sphere(Point3f0(0.0), 1.0f0))
spring   = Gl.Spring()
for i in 1:nspheres
    color = cs[mod1(i, 2)]
    Entity(dio, Gl.Spatial(; position = ps[i]), Gl.Material(), Gl.Shape(),
           Gl.UniformColor(color), Gl.Dynamic(), Gl.Selectable(), sph_geom, spring)
end


text:
function primer33(A)
    A[:, 2:2:end] .= 0
    return A
end

function fakt(n)
    if n == 2
        return n
    end
    return n * fakt(n-1)
end


text:
using Homebrew
using Test

# Print some debugging info
@info("Using Homebrew.jl installed to $(Homebrew.prefix())")

# Restore pkg-config to its installed (or non-installed) state at the end of all of this
pkg_was_installed = Homebrew.installed("pkg-config")
libgfortran_was_installed = Homebrew.installed("staticfloat/juliadeps/libgfortran")

if pkg_was_installed
    @info("Removing pkg-config for our testing...")
    Homebrew.rm("pkg-config")
end

# Add pkg-config
Homebrew.add("pkg-config")
@test Homebrew.installed("pkg-config") == true

# Print versioninfo() to boost coverage
Homebrew.versioninfo()

# Now show that we have it and that it's the right version
function strip_underscores(str)
    range = something(findlast("_", str), 0:0)
    if range.start > 1
        return str[1:range.start-1]
    else
        return str
    end
end
pkgconfig = Homebrew.info("pkg-config")
version = readchomp(`pkg-config --version`)
@test version == strip_underscores(pkgconfig.version)
@test Homebrew.installed(pkgconfig) == true
@info("$(pkgconfig) installed to: $(Homebrew.prefix(pkgconfig))")

@test isdir(Homebrew.prefix("pkg-config"))
@test isdir(Homebrew.prefix(pkgconfig))

# Run through some of the Homebrew API, both with strings and with BrewPkg objects
@test length(filter(x -> x.name == "pkg-config", Homebrew.list())) > 0
@test Homebrew.linked("pkg-config") == true
@test Homebrew.linked(pkgconfig) == true

# Test dependency inspection
@test Homebrew.direct_deps("pkg-config") == []
@test Homebrew.direct_deps(pkgconfig) == []
@test Homebrew.direct_deps("nettle") == [Homebrew.info("gmp")]
@test Homebrew.direct_deps(Homebrew.info("nettle")) == [Homebrew.info("gmp")]

# Run through our sorted deps routines, ensuring that everything is sorted
sortdeps = Homebrew.deps_sorted("pango")
for idx in 1:length(sortdeps)
    for dep in Homebrew.direct_deps(sortdeps[idx])
        depidx = findfirst(x -> (x.name == dep.name), sortdeps)
        @test depidx != 0
        @test depidx < idx
    end
end

# Test that we can probe for bottles properly
@test Homebrew.has_bottle("ack") == false
@test Homebrew.has_bottle("cairo") == true
# I will be a very happy man the day this test starts to fail
@test Homebrew.has_relocatable_bottle("cairo") == false
if Homebrew.has_bottle("staticfloat/juliadeps/libgfortran")
    @test Homebrew.has_relocatable_bottle("staticfloat/juliadeps/libgfortran") == true
end
@test Homebrew.json(pkgconfig)["name"] == "pkg-config"

# Test that has_bottle knows which OSX version we're running on.
@test Homebrew.has_bottle("ld64") == false

# Test that we can translate properly
@info("Translation should pass:")
@test Homebrew.translate_formula("gettext"; verbose=true) == "staticfloat/juliatranslated/gettext"
@info("Translation should fail because it has no bottles:")
@test Homebrew.translate_formula("ack"; verbose=true) == "ack"

if libgfortran_was_installed
    # Remove libgfortran before we start messing around with it
    Homebrew.rm("staticfloat/juliadeps/libgfortran"; force=true)
end

# Make sure translation works properly with other taps
Homebrew.delete_translated_formula("staticfloat/juliadeps/libgfortran"; verbose=true)
@info("Translation should pass because we just deleted libgfortran from translation cache:")
@test Homebrew.translate_formula("staticfloat/juliadeps/libgfortran"; verbose=true) == "staticfloat/juliatranslated/libgfortran"
@info("Translation should fail because libgfortran has already been translated:")
# Do it a second time so we can get coverage of practicing that particular method of bailing out
Homebrew.translate_formula(Homebrew.info("staticfloat/juliadeps/libgfortran"); verbose=true)

# Test that installation of a formula from a tap when it's already been translated works
Homebrew.add("staticfloat/juliadeps/libgfortran"; verbose=true)

if !libgfortran_was_installed
    Homebrew.rm("staticfloat/juliadeps/libgfortran")
end

# Now that we have staticfloat/juliadeps tapped, test to make sure that prefix() works
# with taps properly:
@test Homebrew.prefix("libgfortran") == Homebrew.prefix("staticfloat/juliadeps/libgfortran")

# Test more miscellaneous things
fontconfig = Homebrew.info("staticfloat/juliadeps/fontconfig")
@test Homebrew.formula_path(fontconfig) == joinpath(Homebrew.tappath, "fontconfig.rb")
@test !isempty(Homebrew.read_formula("xz"))
@test !isempty(Homebrew.read_formula(fontconfig))
@info("add() should fail because this actually isn't a package name:")
@test_throws ArgumentError Homebrew.add("thisisntapackagename")

Homebrew.unlink(pkgconfig)
@test Homebrew.installed(pkgconfig) == true
@test Homebrew.linked(pkgconfig) == false
Homebrew.link(pkgconfig)
@test Homebrew.installed(pkgconfig) == true
@test Homebrew.linked(pkgconfig) == true

# Can't really do anything useful with these, but can at least run them to ensure they work
Homebrew.outdated()
Homebrew.update()
Homebrew.postinstall("pkg-config")
Homebrew.postinstall(pkgconfig)
Homebrew.delete_translated_formula("gettext"; verbose=true)
Homebrew.delete_all_translated_formulae(verbose=true)

# Test deletion as well, showing that the array-argument form continues on after errors
Homebrew.rm(pkgconfig)
Homebrew.add(pkgconfig)
@info("rm() should fail because this isn't actually a package name:")
Homebrew.rm(["thisisntapackagename", "pkg-config"])
@test Homebrew.installed("pkg-config") == false
@test Homebrew.linked("pkg-config") == false

if pkg_was_installed
    @info("Adding pkg-config back again...")
    Homebrew.add("pkg-config")
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using JuMP
using MIPLearn
using Cbc

@testset "FileInstance" begin
    @testset "Solve (knapsack)" begin
        data = KnapsackData()
        filename = tempname()
        MIPLearn.save_data(filename, data)
        instance = FileInstance(filename, build_knapsack_model)
        solver = LearningSolver(Cbc.Optimizer)
        solve!(solver, instance)
        h5 = Hdf5Sample("$filename.h5")
        @test h5.get_scalar("mip_wallclock_time") > 0
    end

    @testset "Solve (vpm2)" begin
        data = Dict("filename" => joinpath(@__DIR__, "../fixtures/danoint.mps.gz"))
        build_model(data) = read_from_file(data["filename"])
        filename = tempname()
        MIPLearn.save_data(filename, data)
        instance = FileInstance(filename, build_model)
        solver = LearningSolver(optimizer_with_attributes(Cbc.Optimizer, "seconds" => 1.0))
        solve!(solver, instance)
        h5 = Hdf5Sample("$filename.h5")
        @test h5.get_scalar("mip_wallclock_time") > 0
    end

    @testset "Save and load data" begin
        filename = tempname()
        data = KnapsackData(
            weights = [5.0, 5.0, 5.0],
            prices = [1.0, 1.0, 1.0],
            capacity = 3.0,
        )
        MIPLearn.save_data(filename, data)
        loaded = MIPLearn.load_data(filename)
        @test loaded.weights == [5.0, 5.0, 5.0]
        @test loaded.prices == [1.0, 1.0, 1.0]
        @test loaded.capacity == 3.0
    end
end


text:
"""
    FitnessLayer(
        N::Int64
        )

Layer with `N × N` parameters with dynamics:

```julia
function (L::FitnessLayer)(p)
    f = L.W *p
    ṗ = p .* (f - ones(size(p,1))*p'*f)
    return ( ṗ )
end
```

# Arguments:
- `N::Int64` : layer size on declaration
- `p::Array` : initial condition of dynamics on run
"""
struct FitnessLayer
    W
end

FitnessLayer(N::Int64) = FitnessLayer( zeros(N,N) )
Flux.@functor FitnessLayer

function (L::FitnessLayer)(p)
    f = L.W *p
    ṗ = p .* (f - ones(size(p,1))*p'*f)
    return ( ṗ )
end

"""
    getModel(
        N::Int64
        )

Get compositional Neural ODE (cNODE) for system of size `N`.

# Arguments:
- `N::Int64` : system size
"""
getModel(N::Int64) = NeuralODE(FitnessLayer(N),(0.0,1.0),Tsit5(),saveat=1.0)

"""
    predict(
        cnode::NeuralODE,
        z
        )

Predict composition from normalized collection `z`.

# Arguments:
- `cnode::NeuralODE` : compositional Neural ODE with `FitnessLayer`
- `z::Array` : normalized collection of species with size `N×1`
"""
predict(cnode::NeuralODE,z) = Array(cnode(z).u[end])

_loss(cnode::NeuralODE,z,p) = sum(abs.(p .- predict(cnode,z))) / sum(abs.(p .+ predict(cnode,z)))

"""
    loss(
        node::NeuralODE,
        z,
        p
    )

Compute loss between predictions `q = node(z)` and true data `p`.

# Arguments:
- `node::NeuralODE` : neural ODE with `FitnessLayer`
- `z::Array` : normalized collection of species with size `N×1`
- `p::Array` : true relative abundances, size `N×1`
"""
loss(cnode::NeuralODE,Z,P) = mean([ _loss(cnode,Z[:,i],P[:,i]) for i in 1:size(Z,2) ])

"""
    train_reptile(
        cnode::NeuralODE,
        epochs::Int64,
        mb::Int64,
        LR::Array{Float64},
        Ztrn, Ptrn,
        Zval, Pval,
        Ztst, Ptst,
        report::Int64,
        early_stoping::Int64
    )

Train cNODE using the Reptile meta-learning algorithm.

# Arguments:
- `cnode::NeuralODE` : compositional Neural ODE to train
- `epochs::Int64` : number of epochs for training
- `mb::Int64` : minibatch size
- `LR::Array{Float64}` : array with internal and external learning rates
- `Ztrn, Ptrn` : training data
- `Zval, Pval` : validation data
- `Ztst, Ptst` : test data
- `report::Int64` : interval of epochs to report
- `early_stoping::Int64` : min epochs to start early stopping
"""
function train_reptile(
            cnode::NeuralODE,
            epochs::Int64,
            mb::Int64,
            LR::Array{Float64},
            Ztrn, Ptrn,
            Zval, Pval,
            Ztst, Ptst,
            report::Int64,
            early_stoping::Int64
        )
    loss_train = []
    loss_test = []
    loss_val = []
    stoping = []
    EarlyStoping = []
    es = 0

    W = Flux.params(cnode)
    opt_in = ADAM(LR[1])
    opt_out = ADAM(LR[2])

    for e in 1:epochs
        V = deepcopy(W)
        for (z,p) in eachbatch(shuffleobs((Ztrn,Ptrn)),mb)
            grads = gradient(()->loss(cnode,z,p),W)
            update!(opt_in, W, grads)
        end
        for (w, v) in zip(W,V)
            dv = apply!(opt_out, v, w .- v)
            @. w = v + dv
        end
        push!(loss_train,loss(cnode,Ztrn,Ptrn))
        push!(loss_test,loss(cnode,Ztst,Ptst))
        push!(loss_val,loss(cnode,Zval,Pval))

        push!(stoping,[loss(cnode,Ztrn[:,k],Ptrn[:,k]) for k in 1:size(Ztrn,2)])

        if e > early_stoping
            counts = count(mean(stoping[end-early_stoping:end-1]).<= stoping[end])/size(Ztrn,2)
            push!(EarlyStoping,counts)
            e%report==0 && println(e,
                            "\ttrain - ",round(loss_train[end];digits=3),
                            "\ttest - ",round(loss_test[end];digits=3),
                            '\t',round(counts;digits=3),
                            '\t',round(es;digits=3))
            if counts>0.5
                es +=1
                if es>10
                    "Early stoping at $e..." |> println
                    break
                end
            end
        elseif e>1
            e%report==0 && println(e,
                            "\ttrain - ",round(loss_train[end];digits=3),
                            "\ttest - ",round(loss_test[end];digits=3),'\t',"---")
        end
    end

    return W, loss_train, loss_val, loss_test
end


text:
println("\n\nload_s_earthquake: start")

# number of files to load as sample (0 disables)
# FILE_LIMIT = 2
FILE_LIMIT = 0

using Distributed
if FILE_LIMIT == 0 && nprocs() < length(Sys.cpu_info())
    n = length(Sys.cpu_info()) - nprocs()
    addprocs(n, restrict=true)
end
println("load_s_earthquake: nprocs=$(nprocs())")

@everywhere prepend!(LOAD_PATH, ["Project.toml"])
@everywhere import CSV, LibPQ
@everywhere using Dates, DataFrames, Glob

@everywhere DISASTER_PASSWORD = read("./etc/disaster-pass", String)

FILE_PATHS = glob("*.csv", "./data/usgs_earthquake")
if FILE_LIMIT > 0
    FILE_PATHS = FILE_PATHS[1:FILE_LIMIT]
end
println("load_s_earthquake: FILE_PATHS=$(FILE_PATHS)")

postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

# create empty earthquake staging table
sql =
"""
drop table if exists public.s_earthquake;
create table if not exists public.s_earthquake (
  id                   text primary key,
  time                 text,
  place                text,
  latitude             decimal,
  longitude            decimal,
  depth                decimal,
  magnitude            decimal,
  magnitude_type       text,
  type                 text,
  status               text,
  ogc_fid              integer,
  datestamp            date,
  severity             integer
);

create or replace function public.trigger_s_earthquake_udf()
returns trigger as \$body\$
declare
    _ogc_fid  integer;
begin
    select ogc_fid into _ogc_fid
        from public.tl_2015_us_county
        where st_contains(
            tl_2015_us_county.wkb_geometry,
            ST_SetSRID(ST_MakePoint(new.longitude, new.latitude), 4326));

    -- raise notice 'update_s_earthquake_ogc_fid: id=%, ogc_fid=%', new.id, _ogc_fid;
    new.ogc_fid = _ogc_fid;
    new.datestamp = new.time::date;

    -- use richter scale to create 1-10 severity value
    if new.magnitude < 1.0 then new.severity = 1;
    elsif new.magnitude < 2.0 then new.severity = 2;
    elsif new.magnitude < 3.0 then new.severity = 3;
    elsif new.magnitude < 4.0 then new.severity = 4;
    elsif new.magnitude < 5.0 then new.severity = 5;
    elsif new.magnitude < 6.0 then new.severity = 6;
    elsif new.magnitude < 7.0 then new.severity = 7;
    elsif new.magnitude < 8.0 then new.severity = 8;
    elsif new.magnitude < 9.0 then new.severity = 9;
    else new.severity = 10;
    end if;

    return new;
end; \$body\$ language plpgsql;

drop trigger if exists trigger_s_earthquake on public.s_earthquake;
create trigger trigger_s_earthquake
  before insert on public.s_earthquake
  for each row execute function public.trigger_s_earthquake_udf();
"""
LibPQ.execute(postgres, sql)
LibPQ.close(postgres)

@everywhere function load_s_earthquake(filepath)
    println("\tfilepath=$filepath")
    postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

    df = CSV.read(filepath)
    select!(df, [
        :id, :time, :place, :latitude, :longitude, :depth, :mag, :magType,
        :type, :status, :net, :locationSource, :magSource, :updated
    ])
    select!(df, [
        :id, :time, :place, :latitude, :longitude, :depth, :mag, :magType,
        :type, :status
    ])
    println("\tsize=$(size(df))")
    # println("load_s_earthquake: columns=$(names(df))")

    LibPQ.load!(
        df,
        postgres,
        """
        insert into public.s_earthquake (
            id,
            time,
            place,
            latitude,
            longitude,
            depth,
            magnitude,
            magnitude_type,
            type,
            status
        ) values (
            \$1,
            \$2,
            \$3,
            \$4,
            \$5,
            \$6,
            \$7,
            \$8,
            \$9,
            \$10
        )
        """
    )

    LibPQ.close(postgres)
end

starttime = Dates.now()
jobs = []
for filepath = FILE_PATHS
    job = @spawn load_s_earthquake(filepath)
    push!(jobs, job)
end
for job in jobs
    fetch(job)
end

endtime = Dates.now()
println("load_s_earthquake: elapsed time is $(Dates.canonicalize(Dates.CompoundPeriod(endtime - starttime)))")
println("load_s_earthquake: done")


text:
#=
Utility functions used in the QuantEcon library

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-14

=#

meshgrid(x::Vector, y::Vector) = (repmat(x, 1, length(y))',
                                  repmat(y, 1, length(x)))

# function to return a Range object with points equivalent to calling
# linspace(x_min, x_max, n_x). This is needed because we often use
# CoordInterpGrid from Grid.jl for interpolation and that requires a
# range to be passed as the first argument
linspace_range(x_min, x_max, n_x) = x_min:(x_max - x_min) / (n_x  - 1): x_max


text:
function weighted_mean{T<:Real,W<:Real}(data::Array{T,1}, weights::Array{W,1})
    weights = float(weights) / sum(weights)
    first(float(data)' * weights)
end


text:
## %PREDICTMD_GENERATED_BY%

using PredictMDExtra
PredictMDExtra.import_all()

using PredictMD
PredictMD.import_all()

using CSVFiles
using CategoricalArrays
using DataFrames
using DecisionTree
using Distributions
using FileIO
using GLM
using IterTools
using Knet
using LIBSVM
using LinearAlgebra
using PredictMD
using PredictMDAPI
using PredictMDExtra
using RDatasets
using Random
using StatsModels
using Test
using Unitful

const Schema = StatsModels.Schema

# PREDICTMD IF INCLUDE TEST STATEMENTS
logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
if isnothing(logger)
    logger_stream = devnull
else
    logger_stream = logger.stream
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### Begin project-specific settings

DIRECTORY_CONTAINING_THIS_FILE = @__DIR__
PROJECT_DIRECTORY = dirname(
    joinpath(splitpath(DIRECTORY_CONTAINING_THIS_FILE)...)
    )
PROJECT_OUTPUT_DIRECTORY = joinpath(
    PROJECT_DIRECTORY,
    "output",
    )
mkpath(PROJECT_OUTPUT_DIRECTORY)
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "data"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "models"))
mkpath(joinpath(PROJECT_OUTPUT_DIRECTORY, "plots"))

# PREDICTMD IF INCLUDE TEST STATEMENTS
@debug("PROJECT_OUTPUT_DIRECTORY: ", PROJECT_OUTPUT_DIRECTORY,)
if PredictMD.is_travis_ci()
    PredictMD.cache_to_path!(
        ;
        from = ["cpu_examples", "boston_housing", "output",],
        to = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

### End project-specific settings

### Begin model comparison code

Random.seed!(999)

trainingandtuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_features_df.csv",
    )
trainingandtuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "trainingandtuning_labels_df.csv",
    )
testing_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_features_df.csv",
    )
testing_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "testing_labels_df.csv",
    )
training_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_features_df.csv",
    )
training_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "training_labels_df.csv",
    )
tuning_features_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_features_df.csv",
    )
tuning_labels_df_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "data",
    "tuning_labels_df.csv",
    )
trainingandtuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
trainingandtuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        trainingandtuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )
testing_features_df = DataFrames.DataFrame(
    FileIO.load(
        testing_features_df_filename;
        type_detect_rows = 100,
        )
    )
testing_labels_df = DataFrames.DataFrame(
    FileIO.load(
        testing_labels_df_filename;
        type_detect_rows = 100,
        )
    )
training_features_df = DataFrames.DataFrame(
    FileIO.load(
        training_features_df_filename;
        type_detect_rows = 100,
        )
    )
training_labels_df = DataFrames.DataFrame(
    FileIO.load(
        training_labels_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_features_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_features_df_filename;
        type_detect_rows = 100,
        )
    )
tuning_labels_df = DataFrames.DataFrame(
    FileIO.load(
        tuning_labels_df_filename;
        type_detect_rows = 100,
        )
    )

linear_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "linear_regression.jld2",
    )
random_forest_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "random_forest_regression.jld2",
    )
knet_mlp_regression_filename = joinpath(
    PROJECT_OUTPUT_DIRECTORY,
    "models",
    "knet_mlp_regression.jld2",
    )

# PREDICTMD IF INCLUDE TEST STATEMENTS
linear_regression = nothing
Test.@test isnothing(linear_regression)

random_forest_regression = nothing
Test.@test isnothing(random_forest_regression)

knet_mlp_regression = nothing
Test.@test isnothing(knet_mlp_regression)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

linear_regression =
    PredictMD.load_model(linear_regression_filename)
random_forest_regression =
    PredictMD.load_model(random_forest_regression_filename)
knet_mlp_regression =
    PredictMD.load_model(knet_mlp_regression_filename)

PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)

# PREDICTMD IF INCLUDE TEST STATEMENTS
PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)

PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)

PredictMD.parse_functions!(linear_regression)
PredictMD.parse_functions!(random_forest_regression)
PredictMD.parse_functions!(knet_mlp_regression)
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

all_models = PredictMD.AbstractFittable[
    linear_regression,
    random_forest_regression,
    knet_mlp_regression,
    ]

single_label_name = :MedV

continuous_label_names = Symbol[single_label_name]
categorical_label_names = Symbol[]
label_names = vcat(categorical_label_names, continuous_label_names)

println(logger_stream, "Single label regression metrics, training set: ")
show(
    logger_stream, PredictMD.singlelabelregressionmetrics(
        all_models,
        training_features_df,
        training_labels_df,
        single_label_name,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

println(logger_stream, "Single label regression metrics, testing set: ")
show(
    logger_stream, PredictMD.singlelabelregressionmetrics(
        all_models,
        testing_features_df,
        testing_labels_df,
        single_label_name,
        );
    allrows = true,
    allcols = true,
    splitcols = false,
    )

### End model comparison code

# PREDICTMD IF INCLUDE TEST STATEMENTS
metrics = PredictMD.singlelabelregressionmetrics(all_models,
                                                 testing_features_df,
                                                 testing_labels_df,
                                                 single_label_name)
r2_row = first(
    findall(
        strip.(metrics[:metric]) .== "R^2 (coefficient of determination)"
        )
    )
Test.@test(
    strip(metrics[r2_row, :metric]) == "R^2 (coefficient of determination)"
    )
Test.@test(
    metrics[r2_row, Symbol("Linear regression")] > 0.550
    )
Test.@test(
    metrics[r2_row, Symbol("Random forest")] > 0.550
    )
Test.@test(
    metrics[r2_row, Symbol("Knet MLP")] > 0.300
    )
mse_row = first(
    findall(
        strip.(metrics[:metric]) .== "Mean squared error (MSE)"
        )
    )
Test.@test(
    strip(metrics[mse_row, :metric]) == "Mean squared error (MSE)"
    )
Test.@test(
    metrics[mse_row, Symbol("Linear regression")] < 40.000
    )
Test.@test(
    metrics[mse_row, Symbol("Random forest")] < 40.000
    )
Test.@test(
    metrics[mse_row, Symbol("Knet MLP")] < 65.000
    )
rmse_row = first(
    findall(
        strip.(metrics[:metric]) .== "Root mean square error (RMSE)"
        )
    )
Test.@test(
    strip(metrics[rmse_row, :metric]) == "Root mean square error (RMSE)"
    )
Test.@test(
    metrics[rmse_row, Symbol("Linear regression")] < 6.500
    )
Test.@test(
    metrics[rmse_row, Symbol("Random forest")] < 6.500
    )
Test.@test(
    metrics[rmse_row, Symbol("Knet MLP")] < 8.000
    )
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

# PREDICTMD IF INCLUDE TEST STATEMENTS
if PredictMD.is_travis_ci()
    PredictMD.path_to_cache!(
        ;
        to = ["cpu_examples", "boston_housing", "output",],
        from = [PROJECT_OUTPUT_DIRECTORY],
        )
end
# PREDICTMD ELSE
# PREDICTMD ENDIF INCLUDE TEST STATEMENTS

## %PREDICTMD_GENERATED_BY%


text:
#=
This app creates a simple sidebar layout using inline style arguments and the
dbc_nav component.

dcc_location is used to track the current location, and a callback uses the
current location to render the appropriate page content. The active prop of
each NavLink is set automatically according to the current pathname. To use
this feature you must install dash-bootstrap-components >= 0.11.0.

For more details on building multi-page Dash applications, check out the Dash
documentation: https://dash.plot.ly/urls
=#

using Dash, DashBootstrapComponents, DashCoreComponents, DashHtmlComponents

app = dash(external_stylesheets = [dbc_themes.BOOTSTRAP])

# the style arguments for the sidebar. We use position => fixed and a fixed width
SIDEBAR_STYLE = Dict(
    "position" => "fixed",
    "top" => 0,
    "left" => 0,
    "bottom" => 0,
    "width" => "16rem",
    "padding" => "2rem 1rem",
    "background-color" => "#f8f9fa",
)

# the styles for the main content position it to the right of the sidebar and
# add some padding.
CONTENT_STYLE =
    Dict("margin-left" => "18rem", "margin-right" => "2rem", "padding" => "2rem 1rem")

sidebar = html_div(
    [
        html_h2("Sidebar", className = "display-4"),
        html_hr(),
        html_p("A simple sidebar layout with navigation links", className = "lead"),
        dbc_nav(
            [
                dbc_navlink("Home", href = "/", active = "exact"),
                dbc_navlink("Page 1", href = "/page-1", active = "exact"),
                dbc_navlink("Page 2", href = "/page-2", active = "exact"),
            ],
            vertical = true,
            pills = true,
        ),
    ],
    style = SIDEBAR_STYLE,
)

content = html_div(id = "page-content", style = CONTENT_STYLE)

app.layout = html_div([dcc_location(id = "url"), sidebar, content])

callback!(app, Output("page-content", "children"), Input("url", "pathname")) do pathname
    if pathname == "/"
        return html_p("This is the content of the home page!")
    elseif pathname == "/page-1"
        return html_p("This is the content of page 1. Yay!")
    elseif pathname == "/page-2"
        return html_p("Oh cool, this is page 2!")
    else
        # If the user tries to reach a different page, return a 404 message
        return dbc_jumbotron([
            html_h1("404: Not found", className = "text-danger"),
            html_hr(),
            html_p("The pathname $pathname was not recognised..."),
        ])
    end
end

run_server(app, "0.0.0.0", 8888)

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components