text:
################################################################################
# FourIntersectionRoadwayOptions
################################################################################

@with_kw mutable struct FourIntersectionRoadwayOptions{T} <: RoadwayOptions
    # Options
	"Length of the lanes."
	lane_length::T=4.0

	"Width of the lanes."
	lane_width::T=0.4

	"Turn radius."
	turn_radius::T=0.1
end

function build_roadway(opts::FourIntersectionRoadwayOptions{T}) where {T}
	ll = opts.lane_length
	lw = opts.lane_width
	tr = opts.turn_radius
	#                  x5      x17       x6
	#                  |        .        |
	#                  |        .        |
	#                  |        .        |
	#                  |        .        |
	#            x29   x25      .        x26  x30
	#                  |        .        |
	# x1_________x33___|x9     x18      x10___x22___________x2
	#
	#                         (0,0)
	# x13..............x14	   x21      x15.................x16
	#
	#
	# x3_________x23___x11     x19      x12___x24___________x4
	#                  |        .        |
	#            x31   x27        .      x28  x32
	#                  |        .        |
	#                  |        .        |
	#                  |        .        |
	#                  |        .        |
	#                  |        .        |
	#                  x7      x20       x8
	#

	x1 = [ -ll/2, lw/2]
	x2 = [ ll/2, lw/2]
	x3 = [ -ll/2, -lw/2]
	x4 = [ ll/2, -lw/2]
	x5 = [ -lw/2, ll/2]
	x6 = [ lw/2, ll/2]
	x7 = [ -lw/2, -ll/2]
	x8 = [ lw/2, -ll/2]
	x9 = [ -lw/2, lw/2]
	x10 = [ lw/2, lw/2]
	x11 = [ -lw/2, -lw/2]
	x12 = [ lw/2, -lw/2]
	x13 = [ -ll/2, 0.]
	x14 = [ -lw/2, 0.]
	x15 = [ lw/2, 0.]
	x16 = [ ll/2, 0.]
	x17 = [ 0., ll/2]
	x18 = [ 0., lw/2]
	x19 = [ 0., -lw/2]
	x20 = [ 0., -ll/2]

	x21 = [0., 0.]
	x22 = [+lw/2 + tr, +lw/2]
	x23 = [-lw/2 - tr, -lw/2]
	x24 = [+lw/2 + tr, -lw/2]
	x25 = [-lw/2, +lw/2 + tr]
	x26 = [+lw/2, +lw/2 + tr]
	x27 = [-lw/2, -lw/2 - tr]
	x28 = [+lw/2, -lw/2 - tr]
	x29 = [-lw/2 - tr, +lw/2 + tr]
	x30 = [+lw/2 + tr, +lw/2 + tr]
	x31 = [-lw/2 - tr, -lw/2 - tr]
	x32 = [+lw/2 + tr, -lw/2 - tr]
	x33 = [-lw/2 - tr, +lw/2]

	v = [0., 1.]
	w = [1., 0.]
    w1 = Wall(x1, x9, v)
	w2 = Wall(x10, x2, v)
	w3 = Wall(x3, x11, -v)
	w4 = Wall(x12, x4, -v)
	w5 = Wall(x9, x5, -w)
	w6 = Wall(x7, x11, -w)
	w7 = Wall(x8, x12,  w)
	w8 = Wall(x10, x6,  w)
	w9 = Wall(x13, x14, v)
	w10 = Wall(x13, x14, -v)
	w11 = Wall(x15, x16,  v)
	w12 = Wall(x15, x16, -v)
	w13 = Wall(x18, x17, -w)
	w14 = Wall(x18, x17,  w)
	w15 = Wall(x20, x19, -w)
	w16 = Wall(x20, x19,  w)

	# West
	trv = [tr, 0]
	trw = [0, tr]
	w101  = Wall(x3,  x4,  -v)
	w102  = Wall(x13, x16,  v)

	w103  = Wall(x3,  x23, -v)
	w104  = Wall(x13, x21,  v)
	w105  = Wall(x7,  x27, -w)
	w106  = Wall(x20, x21,  w)
	c107  = CircularWall(x31..., tr)

	w108  = Wall(x3,  x12, -v)
	w109  = Wall(x13, x21-trv,  v)
	w110  = Wall(x17, x21+trw, -w)
	w111  = Wall(x6 , x12,  w)
	c112  = CircularWall((x21-trw+trv)..., tr)

	# East
	w201  = Wall(x1,  x2,  v)
	w202  = Wall(x13, x16,  -v)

	w203  = Wall(x22,  x2, v)
	w204  = Wall(x16, x21,  -v)
	w205  = Wall(x26,  x6, w)
	w206  = Wall(x21, x17,  -w)
	c207  = CircularWall(x30..., tr)

	w208  = Wall(x9,  x2, v)
	w209  = Wall(x21+trv, x16, -v)
	w210  = Wall(x21-trw, x20,  w)
	w211  = Wall(x7 , x9,  -w)
	c212  = CircularWall((x21+trw-trv)..., tr)

	# North
	w301  = Wall(x5,  x7,  -w)
	w302  = Wall(x17, x20,  w)

	w303  = Wall(x25,  x5, -w)
	w304  = Wall(x17, x21,  w)
	w305  = Wall(x1,  x33, v)
	w306  = Wall(x21, x13,  -v)
	c307  = CircularWall(x29..., tr)

	w308  = Wall(x5,  x11, -w)
	w309  = Wall(x21+trw, x17, w)
	w310  = Wall(x21+trv, x16,  v)
	w311  = Wall(x4, x11,  -v)
	c312  = CircularWall((x21-trw-trv)..., tr)

	# South
	w401  = Wall(x6,  x8,  w)
	w402  = Wall(x17, x20,  -w)

	w403  = Wall(x8,  x28, w)
	w404  = Wall(x20, x21,  -w)
	w405  = Wall(x4,  x24, -v)
	w406  = Wall(x21, x16,  v)
	c407  = CircularWall(x32..., tr)

	w408  = Wall(x8,  x10,  w)
	w409  = Wall(x21-trw, x20, -w)
	w410  = Wall(x21-trv, x13,  -v)
	w411  = Wall(x1, x10,  v)
	c412  = CircularWall((x21+trw+trv)..., tr)

	f_west_straight(s,v=0.) = VehicleState(s-ll/2, -lw/4, 0., v)
	function f_west_right(s,v=0.)
		if s <= ll/2-lw/4
			vs = VehicleState(s-ll/2, -lw/4, 0., v)
		else
			vs = VehicleState(-lw/4, -lw/4-(s-(ll/2-lw/4)), -π/2, v)
		end
		return vs
	end
	function f_west_left(s,v=0.)
		if s <= ll/2+lw/4
			vs = VehicleState(s-ll/2, -lw/4, 0., v)
		else
			vs = VehicleState(+lw/4, -lw/4+(s-(ll/2+lw/4)), π/2, v)
		end
		return vs
	end

	f_east_straight(s,v=0.) = VehicleState(ll/2 - s, lw/4, 3.14, v)
	function f_east_right(s,v=0.)
		if s <= ll/2-lw/4
			vs = VehicleState(ll/2 - s, lw/4, 3.14, v)
		else
			vs = VehicleState(lw/4, lw/4+(s-(ll/2-lw/4)), π/2, v)
		end
		return vs
	end
	function f_east_left(s,v=0.)
		if s <= ll/2 + lw/4
			vs = VehicleState(ll/2 - s, lw/4, 3.14, v)
		else
			vs = VehicleState(-lw/4, lw/4-(s-(ll/2+lw/4)), -1.57, v)
		end
		return vs
	end
	f_north_straight(s,v=0.) = VehicleState(-lw/4, ll/2 - s, -1.57, v)
	function f_north_right(s,v=0.)
		if s <= ll/2-lw/4
			vs = VehicleState(-lw/4, ll/2 - s, -1.57, v)
		else
			vs = VehicleState(-lw/4, lw/4+(s-(ll/2-lw/4)), -π, v)
		end
		return vs
	end
	function f_north_left(s,v=0.)
		if s <= ll/2 + lw/4
			vs = VehicleState(-lw/4, ll/2 - s, -1.57, v)
		else
			vs = VehicleState(-lw/4, -lw/4+(s-(ll/2+lw/4)), -π, v)
		end
		return vs
	end
	f_south_straight(s,v=0.) = VehicleState(lw/4, s - ll/2, 1.57, v)
	function f_south_right(s,v=0.)
		if s <= ll/2-lw/4
			vs = VehicleState(lw/4, s - ll/2, 1.57, v)
		else
			vs = VehicleState(lw/4, -lw/4+(s-(ll/2-lw/4)), 0., v)
		end
		return vs
	end
	function f_south_left(s,v=0.)
		if s <= ll/2 + lw/4
			vs = VehicleState(lw/4, s - ll/2, 1.57, v)
		else
			vs = VehicleState(lw/4, +lw/4+(s-(ll/2+lw/4)), 3.14, v)
		end
		return vs
	end

	# <start_direction>_lane
	# West lane straight
	west_lane_straight = Lane(1,
					:west_lane,
					[w101, w102],
					Vector{CircularWall}(),
					StartingArea(
						VehicleState(-ll/2, -lw/4, 0.0,  0.00),
						VehicleState(-ll/2, -lw/4, 0.0, -0.05),
						VehicleState(-ll/2, -lw/4, 0.4,  0.05),
						),
					f_west_straight,
					)
	# West lane with right turn
	west_lane_right = Lane(2,
					:west_lane_right,
					[w103, w104, w105, w106],
					[c107],
					StartingArea(
						VehicleState(-ll/2, -lw/4, 0.0,  0.00),
						VehicleState(-ll/2, -lw/4, 0.0, -0.05),
						VehicleState(-ll/2, -lw/4, 0.4,  0.05),
						),
					f_west_right,
					)
	# West lane with left turn
	west_lane_left = Lane(3,
					:west_lane_left,
					[w108, w109, w110, w111],
					[c112],
					StartingArea(
						VehicleState(-ll/2, -lw/4, 0.0,  0.00),
						VehicleState(-ll/2, -lw/4, 0.0, -0.05),
						VehicleState(-ll/2, -lw/4, 0.4,  0.05),
						),
					f_west_left,
					)

	# East line straight
	east_lane_straight = Lane(4,
					:east_lane,
					[w201, w202],
					Vector{CircularWall}(),
					StartingArea(
						VehicleState(ll/2, lw/4, 3.14,  0.00),
						VehicleState(ll/2, lw/4, 3.10, -0.05),
						VehicleState(ll/2, lw/4, 3.18,  0.05),
						),
					f_east_straight,
					)
	# East lane with right turn
	east_lane_right = Lane(5,
					:east_lane_right,
					[w203, w204, w205, w206],
					[c207],
					StartingArea(
						VehicleState(ll/2, lw/4, 3.14,  0.00),
						VehicleState(ll/2, lw/4, 3.10, -0.05),
						VehicleState(ll/2, lw/4, 3.18,  0.05),
						),
					f_east_right,
					)
	# East lane with left turn
	east_lane_left = Lane(6,
					:east_lane_left,
					[w208, w209, w210, w211],
					[c212],
					StartingArea(
						VehicleState(ll/2, lw/4, 3.14,  0.00),
						VehicleState(ll/2, lw/4, 3.10, -0.05),
						VehicleState(ll/2, lw/4, 3.18,  0.05),
						),
					f_east_left,
					)
	# North lane straigth
	north_lane_straight = Lane(7,
					:north_lane,
					[w301, w302],
					Vector{CircularWall}(),
					StartingArea(
						VehicleState(-lw/4, ll/2, -1.57,  0.00),
						VehicleState(-lw/4, ll/2, -1.17, -0.05),
						VehicleState(-lw/4, ll/2, -1.97,  0.05),
						),
					f_north_straight,
					)
	# North lane with right turn
	north_lane_right = Lane(8,
					:north_lane_right,
					[w303, w304, w305, w306],
					[c307],
					StartingArea(
						VehicleState(-lw/4, ll/2, -1.57,  0.00),
						VehicleState(-lw/4, ll/2, -1.17, -0.05),
						VehicleState(-lw/4, ll/2, -1.97,  0.05),
						),
					f_north_right,
					)
	# North lane with left turn
	north_lane_left = Lane(9,
					:north_lane_left,
					[w308, w309, w310, w311],
					[c312],
					StartingArea(
						VehicleState(-lw/4, ll/2, -1.57,  0.00),
						VehicleState(-lw/4, ll/2, -1.17, -0.05),
						VehicleState(-lw/4, ll/2, -1.97,  0.05),
						),
					f_north_left,
					)
	# South lane straigth
	south_lane_straight = Lane(10,
					:south_lane,
					[w401, w402],
					Vector{CircularWall}(),
					StartingArea(
						VehicleState(lw/4, -ll/2, 1.57,  0.00),
						VehicleState(lw/4, -ll/2, 1.17, -0.05),
						VehicleState(lw/4, -ll/2, 1.97,  0.05),
						),
					f_south_straight,
					)
	# North lane with right turn
	south_lane_right = Lane(11,
					:north_lane_right,
					[w403, w404, w405, w406],
					[c407],
					StartingArea(
						VehicleState(lw/4, -ll/2, 1.57,  0.00),
						VehicleState(lw/4, -ll/2, 1.17, -0.05),
						VehicleState(lw/4, -ll/2, 1.97,  0.05),
						),
					f_north_right,
					)
	# North lane with left turn
	south_lane_left = Lane(12,
					:south_lane_left,
					[w408, w409, w410, w411],
					[c412],
					StartingArea(
						VehicleState(lw/4, -ll/2, 1.57,  0.00),
						VehicleState(lw/4, -ll/2, 1.17, -0.05),
						VehicleState(lw/4, -ll/2, 1.97,  0.05),
						),
					f_south_left,
					)

	lanes = [
		west_lane_straight, west_lane_right, west_lane_left,
		east_lane_straight, east_lane_right, east_lane_left,
		north_lane_straight, north_lane_right, north_lane_left,
		south_lane_straight, south_lane_right, south_lane_left]
	roadway = Roadway(lanes, opts)
	return roadway
end


text:
typealias Shape Matrix{Float64}

type ModelFrame
    mini::Int
    minj::Int
    maxi::Int
    maxj::Int
    h::Int
    w::Int
    ModelFrame(mini::Int, minj::Int, maxi::Int, maxj::Int) =
        new(mini, minj, maxi, maxj, @compat(Int(ceil(maxi - mini + 3))), @compat(Int(ceil(maxj - minj + 3))))
end

Base.size(frame::ModelFrame) = (frame.h, frame.w)

type Grad2D
    di::Matrix{Float64}
    dj::Matrix{Float64}
end

Base.show(io::IO, grad::Grad2D) = print(io, "Grad2D$(size(grad.di))")

type AAModel
    np::Int
    nc::Int
    frame::ModelFrame
    wparams::PAWarpParams
    s0::Vector{Float64}
    s_star::Matrix{Float64}
    S::Matrix{Float64}
    A0::Vector{Float64}
    A::Matrix{Float64}
    dA0::Vector{Grad2D}
    dW_dp::Array{Float64, 4}
    dN_dq::Array{Float64, 4}
    SD::Matrix{Float64}
    H::Matrix{Float64}
    invH::Matrix{Float64}
    R::Matrix{Float64}
    AAModel() = new()
end

Base.show(io::IO, m::AAModel) = print(io, "AAModel($(m.np))")

ext: jl


text:
@testset "Homotopies" begin
    @testset "ParameterHomotopy" begin
        @var x y a b c

        f = [(2 * x^2 + b^2 * y^3 + 2 * a * x * y)^3, (a + c)^4 * x + y^2]
        F = System(f, [x, y], [a, b, c])

        p = [5.2, -1.3, 9.3]
        q = [2.6, 3.3, 2.3]

        H = HC.ParameterHomotopy(F, p, q)

        v = [0.192, 2.21]
        t = 0.232

        u = zeros(ComplexF64, 2)
        U = zeros(ComplexF64, 2, 2)

        HC.evaluate!(u, H, v, t)
        @test u ≈ f([x, y] => v, [a, b, c] => t * p + (1 - t) * q)

        HC.taylor!(u, Val(1), H, TaylorVector{1}(Matrix(v')), t)
        @test u ≈ let
            @var s sp[1:3] sq[1:3]
            pt = s .* sp .+ (1 .- s) .* sq
            differentiate(subs(f, [a, b, c] => pt), s)(
                [x, y] => v,
                sp => p,
                sq => q,
                s => t,
            )
        end

        HC.evaluate_and_jacobian!(u, U, H, v, t)
        @test U ≈ differentiate(f, [x, y])([x, y] => v, [a, b, c] => t * p + (1 - t) * q)
    end

    @testset "ModelKitHomotopy" begin
        @var x y a b c
        f = [(2 * x^2 + b^2 * y^3 + 2 * a * x * y)^3, (a + c)^4 * x + y^2]

        @var s sp[1:3] sq[1:3]
        g = subs(f, [a, b, c] => s .* sp .+ (1 .- s) .* sq)

        H = Homotopy(g, [x, y], s, [sp; sq])

        p = [5.2, -1.3, 9.3]
        q = [2.6, 3.3, 2.3]

        H = HC.ModelKitHomotopy(H, [p; q])
        @test size(H) == (2, 2)

        v = [0.192, 2.21]
        t = 0.232

        u = zeros(ComplexF64, 2)
        U = zeros(ComplexF64, 2, 2)

        HC.evaluate!(u, H, v, t)
        @test u ≈ f([x, y] => v, [a, b, c] => t * p + (1 - t) * q)

        HC.taylor!(u, Val(1), H, TaylorVector{1}(Matrix(v')), t)
        @test u ≈ let
            @var s sp[1:3] sq[1:3]
            pt = s .* sp .+ (1 .- s) .* sq
            differentiate(subs(f, [a, b, c] => pt), s)(
                [x, y] => v,
                sp => p,
                sq => q,
                s => t,
            )
        end

        HC.evaluate_and_jacobian!(u, U, H, v, t)
        @test U ≈ differentiate(f, [x, y])([x, y] => v, [a, b, c] => t * p + (1 - t) * q)

        H = ModelKit.Homotopy([(2 * x^2 + y^3 + 2 * a * y)^3, x + y^2], [x, y], a)
        @test ModelKitHomotopy(ModelKit.compile(H)) isa ModelKitHomotopy
    end
end


text:
# SINGLE NODE/CORE VERSION
# Filtering out less than in ltw, since this corpus small

import Pkg; Pkg.instantiate()

include("/users/yh31/scratch/projects/gigaword_64k/src/gigaword_64k.jl")

const path_ltw_dir = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga5/data/ltw_eng"

const path_intermed_data = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga_topk/intermediate_files"
const path_intermed_ltw = joinpath(path_intermed_data, "ltw")

# const path_of_tree = path_ltw_dir
# const path_output = path_intermed_ltw
# const n_items = 8

process_part_of_tree(path_ltw_dir, path_intermed_ltw, 6)
# Sat 6pm: OK, calling the process_part_of_tree function didn't work, even though I had used an include and one less worker than i'd requested from slurm. one of the workers became uncommunicative for some reason

# So let's try putting the filetrees stuff in here...

#=
data_tree = FileTree(path_of_tree)

# 1. Get the list of names of files in tree
filenames = name.(files(data_tree))

# 2. For each year that's represented, randomly choose min(array_size, n_items) of the files in that year

pick_elts_from_arr(arr) = pick_elts(arr, n_items)

chosen_files = @_ filenames |>
       collect(IterTools.groupby(year_from_fnm, __)) |>
       map(pick_elts_from_arr, __) |>
       vcat(__...) |> Set(__)

# 3. Get the resulting tree with `filter`
filtered_tree = filter(x->x.name ∈ chosen_files, data_tree, dirs=false)

# 4. load it; count words; save it
loaded_tree = FileTrees.load(filtered_tree; lazy = true) do file
                  @_ file |>
                  string(FileTrees.path(__)) |>
                  read_and_wc(__) |>
                  df_from_acc(__)
             end

out_tree = FileTrees.rename(loaded_tree, path_output)

FileTrees.save(out_tree) do acc_df
   savejdf(string(FileTrees.path(acc_df)) * ".jdf", FileTrees.get(acc_df))
end

#       @warn "failed to save $(string(FileTrees.path(acc_df)))"
#   end
=#


text:
using ReverseSearch

function reverse_search_params_for_complete_graph(n)
    function adj(vertex, index)
        if n != vertex
            index
        else
            Union{}
        end
    end

    max_degree = n
    sink = 1

    function f(vertex)
        if vertex == 1
            Union{}
        else
            vertex - 1
        end
    end

    ReverseSearchParams(adj, max_degree, sink, f)
end

function reverse_search_vertices_for_complete_graph(n)
    reverse_search_params = reverse_search_params_for_complete_graph(n)
    reverse_search_task = @task reverse_search_producer(reverse_search_params)
    sort(collect(reverse_search_task))
end


text:
# Example of vOptGeneric, for a 2IP (the 2ap03.mop is a 2LAP)
# the numerical instance is read on a file, according to the MOP format
# July 2017-updated July 2021

# ---- Packages to use
using vOptGeneric, JuMP, GLPK

m = parseMOP( "2ap03.mop" , GLPK.Optimizer ) # comment: set up the path to the file if required

# ---- Invoking the solver (Chalmet method)
@time vSolve( m, method=:chalmet, step=0.5 )

printX_E(m)
getY_N(m)

max_stars_repo_name: vOptSolver/vOptGeneric

max_issues_repo_name: vOptSolver/vOptGeneric

max_forks_repo_name: vOptSolver/vOptGeneric


text:
function autotune(nlp_raw, pars::Class_parameters)
    nlp = Class_CUTEst(nlp_raw);

    best_it = Inf
    best_mu_ratio = 0.0
    success = true

    for i = -2:5
      local_timer = class_advanced_timer()
      start_advanced_timer(local_timer)
      mu_ratio = 10.0^(i)
      pars.init.mu_scale = mu_ratio
      println("--------  MU_SCALE = $(mu_ratio) ----------")
      start_advanced_timer(local_timer, "INIT");
      intial_it = init(nlp, pars, local_timer);
      pause_advanced_timer(local_timer, "INIT");

      @assert(is_feasible(intial_it, pars.ls.comp_feas))
      iter, status, hist, t, err = one_phase_IPM(intial_it, pars, local_timer);
      if t < best_it - 1 && status == :Optimal
        best_it = t
        best_mu_ratio = mu_ratio
      end

      pause_advanced_timer(local_timer)
      print_timer_stats(local_timer)
    end

    return success, best_it, best_mu_ratio
end


text:
# MadNLP.jl
# Modified from Ipopt.jl (https://github.com/jump-dev/Ipopt.jl)

@kwdef mutable struct VariableInfo
    lower_bound::Float64 = -Inf
    has_lower_bound::Bool = false
    lower_bound_dual_start::Union{Nothing, Float64} = nothing
    upper_bound::Float64 = Inf
    has_upper_bound::Bool = false
    upper_bound_dual_start::Union{Nothing, Float64} = nothing
    is_fixed::Bool = false
    start::Union{Nothing,Float64} = nothing
end

mutable struct ConstraintInfo{F, S}
    func::F
    set::S
    dual_start::Union{Nothing, Float64}
end
ConstraintInfo(func, set) = ConstraintInfo(func, set, nothing)

mutable struct Optimizer <: MOI.AbstractOptimizer
    ips::Union{Nothing,InteriorPointSolver}
    nlp::Union{Nothing,AbstractNLPModel}
    result::Union{Nothing,MadNLPExecutionStats{Float64}}

    # Problem data.
    variable_info::Vector{VariableInfo}
    nlp_data::MOI.NLPBlockData
    sense::MOI.OptimizationSense
    objective::Union{
        MOI.VariableIndex,MOI.ScalarAffineFunction{Float64},MOI.ScalarQuadraticFunction{Float64},Nothing}

    linear_le_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}}
    linear_ge_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}}}
    linear_eq_constraints::Vector{ConstraintInfo{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}}}
    quadratic_le_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64}}}
    quadratic_ge_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64},MOI.GreaterThan{Float64}}}
    quadratic_eq_constraints::Vector{ConstraintInfo{MOI.ScalarQuadraticFunction{Float64},MOI.EqualTo{Float64}}}
    nlp_dual_start::Union{Nothing,Vector{Float64}}

    # Parameters.
    option_dict::Dict{Symbol, Any}

    # Solution attributes.
    solve_time::Float64
end

function Optimizer(;kwargs...)
    option_dict = Dict{Symbol, Any}()
    for (name, value) in kwargs
        option_dict[name] = value
    end
    return Optimizer(nothing,nothing,nothing,[],empty_nlp_data(),MOI.FEASIBILITY_SENSE,
                     nothing,[],[],[],[],[],[],nothing,option_dict,NaN)
end

# define empty nlp evaluator
struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
MOI.eval_constraint(::EmptyNLPEvaluator, g, x) = @assert length(g) == 0
MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x) = fill!(g, 0.0)
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int,Int}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int,Int}[]
MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x) = nothing
MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, σ, μ) = nothing
empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

# for throw_if_valid
MOI.is_valid(model::Optimizer, vi::MOI.VariableIndex) = vi.value in eachindex(model.variable_info)
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}})
    vi = MOI.VariableIndex(ci.value)
    return MOI.is_valid(model, vi) && has_upper_bound(model, vi)
end
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}})
    vi = MOI.VariableIndex(ci.value)
    return MOI.is_valid(model, vi) && has_lower_bound(model, vi)
end
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}})
    vi = MOI.VariableIndex(ci.value)
    return MOI.is_valid(model, vi) && is_fixed(model, vi)
end

# supported obj/var/cons
MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
MOI.supports(::Optimizer,::MOI.ObjectiveFunction{MOI.VariableIndex}) = true
MOI.supports(::Optimizer,::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supports(::Optimizer,::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,::Type{MOI.VariableIndex}) = true
MOI.supports_constraint(::Optimizer,::Type{MOI.VariableIndex},::Type{MOI.LessThan{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.VariableIndex},::Type{MOI.GreaterThan{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.VariableIndex},::Type{MOI.EqualTo{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.LessThan{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.GreaterThan{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.EqualTo{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarQuadraticFunction{Float64}},::Type{MOI.LessThan{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarQuadraticFunction{Float64}},::Type{MOI.GreaterThan{Float64}})=true
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarQuadraticFunction{Float64}},::Type{MOI.EqualTo{Float64}})=true

## TODO
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarAffineFunction{Float64}},::Type{MOI.Interval{Float64}})=false
MOI.supports_constraint(::Optimizer,::Type{MOI.ScalarQuadraticFunction{Float64}},::Type{MOI.Interval{Float64}})=false

MOI.get(::Optimizer, ::MOI.SolverName) = "MadNLP"
MOI.get(::Optimizer, ::MOI.SolverVersion) = version()
MOI.get(model::Optimizer,::MOI.ObjectiveFunctionType)=typeof(model.objective)
MOI.get(model::Optimizer,::MOI.NumberOfVariables)=length(model.variable_info)
MOI.get(model::Optimizer,::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}})=length(model.linear_le_constraints)
MOI.get(model::Optimizer,::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}})=length(model.linear_eq_constraints)
MOI.get(model::Optimizer,::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}})=length(model.linear_ge_constraints)
MOI.get(model::Optimizer,::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.LessThan{Float64}})=count(e->e.has_upper_bound,model.variable_info)
MOI.get(model::Optimizer,::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.EqualTo{Float64}})=count(e->e.is_fixed,model.variable_info)
MOI.get(model::Optimizer,::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.GreaterThan{Float64}})=count(e->e.has_lower_bound,model.variable_info)
MOI.get(model::Optimizer,::MOI.ObjectiveFunction) = model.objective
MOI.get(model::Optimizer,::MOI.ListOfVariableIndices) = [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
MOI.get(model::Optimizer,::MOI.BarrierIterations) = model.ips.cnt.k

macro define_get_con_attr(attr,function_type, set_type, prefix, attrfield)
    array_name = Symbol(string(prefix) * "_constraints")
    quote
        function MOI.get(model::Optimizer, ::$attr,
                         c::MOI.ConstraintIndex{$function_type,$set_type})
            return model.$(array_name)[c.value].$attrfield
        end
    end
end

@define_get_con_attr(MOI.ConstraintSet,MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64},linear_le,set)
@define_get_con_attr(MOI.ConstraintSet,MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64},linear_ge,set)
@define_get_con_attr(MOI.ConstraintSet,MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64},linear_eq,set)
@define_get_con_attr(MOI.ConstraintSet,MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64},quadratic_le,set)
@define_get_con_attr(MOI.ConstraintSet,MOI.ScalarQuadraticFunction{Float64},MOI.GreaterThan{Float64},quadratic_ge,set)
@define_get_con_attr(MOI.ConstraintSet,MOI.ScalarQuadraticFunction{Float64},MOI.EqualTo{Float64},quadratic_eq,set)

@define_get_con_attr(MOI.ConstraintFunction,MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64},linear_le,func)
@define_get_con_attr(MOI.ConstraintFunction,MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64},linear_ge,func)
@define_get_con_attr(MOI.ConstraintFunction,MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64},linear_eq,func)

function MOI.get(model::Optimizer,::MOI.ConstraintSet,
                 c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}})
    return MOI.LessThan{Float64}(model.variable_info[c.value].upper_bound)
end

function MOI.get(model::Optimizer,::MOI.ConstraintSet,
                 c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}})
    return MOI.EqualTo{Float64}(model.variable_info[c.value].lower_bound)
end

function MOI.get(model::Optimizer,::MOI.ConstraintSet,
                 c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}})
    return MOI.GreaterThan{Float64}(model.variable_info[c.value].lower_bound)
end
function MOI.get(model::Optimizer,::MOI.ConstraintFunction,
                 c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}})
    return MOI.VariableIndex(c.value)
end
function MOI.get(model::Optimizer,::MOI.ConstraintFunction,
                 c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}})
    return MOI.VariableIndex(c.value)
end
function MOI.get(
    model::Optimizer,::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}})
    return MOI.VariableIndex(c.value)
end

function MOI.get(model::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    constraints = Set{Tuple{DataType, DataType}}()

    for info in model.variable_info
        info.has_lower_bound && push!(constraints, (MOI.VariableIndex, MOI.LessThan{Float64}))
        info.has_upper_bound && push!(constraints, (MOI.VariableIndex, MOI.GreaterThan{Float64}))
        info.is_fixed && push!(constraints, (MOI.VariableIndex, MOI.EqualTo{Float64}))
    end

    isempty(model.linear_le_constraints) ||
        push!(constraints,(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}))
    isempty(model.linear_ge_constraints) ||
        push!(constraints, (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}))
    isempty(model.linear_eq_constraints) ||
        push!(constraints, (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}))
    isempty(model.quadratic_le_constraints) ||
        push!(constraints, (MOI.ScalarQuadraticFunction{Float64}, MOI.LessThan{Float64}))
    isempty(model.quadratic_ge_constraints) ||
        push!(constraints, (MOI.ScalarQuadraticFunction{Float64}, MOI.GreaterThan{Float64}))
    isempty(model.quadratic_eq_constraints) ||
        push!(constraints, (MOI.ScalarQuadraticFunction{Float64}, MOI.EqualTo{Float64}))

    return collect(constraints)
end
function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}
)
    return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}.(eachindex(model.linear_le_constraints))
end
function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
)
    return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}.(eachindex(model.linear_eq_constraints))
end
function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}
)
    return MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}.(eachindex(model.linear_ge_constraints))
end
function MOI.get(model::Optimizer, ::MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.LessThan{Float64}})
    dict = Dict(model.variable_info[i] => i for i in 1:length(model.variable_info))
    filter!(info -> info.first.has_upper_bound, dict)
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}.(values(dict))
end
function MOI.get(model::Optimizer, ::MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.EqualTo{Float64}})
    dict = Dict(model.variable_info[i] => i for i in 1:length(model.variable_info))
    filter!(info -> info.first.is_fixed, dict)
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}.(values(dict))
end
function MOI.get(model::Optimizer, ::MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.GreaterThan{Float64}})
    dict = Dict(model.variable_info[i] => i for i in 1:length(model.variable_info))
    filter!(info -> info.first.has_lower_bound, dict)
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}.(values(dict))
end

function MOI.set(model::Optimizer, ::MOI.ConstraintSet,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}},
                 set::MOI.LessThan{Float64})
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].upper_bound = set.upper
    return
end

function MOI.set(model::Optimizer, ::MOI.ConstraintSet,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}},
                 set::MOI.GreaterThan{Float64})
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = set.lower
    return
end

function MOI.set(model::Optimizer, ::MOI.ConstraintSet,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}},
                 set::MOI.EqualTo{Float64})
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = set.value
    model.variable_info[ci.value].upper_bound = set.value
    return
end

# default_copy_to
MOI.supports_incremental_interface(model::Optimizer) = true
MOI.copy_to(model::Optimizer, src::MOI.ModelLike) = MOIU.default_copy_to(model, src)

# objective sense
MOI.supports(::Optimizer,::MOI.ObjectiveSense) = true
MOI.set(model::Optimizer, ::MOI.ObjectiveSense,sense::MOI.OptimizationSense) = (model.sense = sense)
MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.sense

# silence
const SILENT_KEY = :print_level
const SILENT_VAL = ERROR
MOI.supports(::Optimizer,::MOI.Silent) = true
MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)= value ?
    MOI.set(model,MOI.RawOptimizerAttribute(String(SILENT_KEY)),SILENT_VAL) : delete!(model.option_dict,SILENT_KEY)
MOI.get(model::Optimizer, ::MOI.Silent) = get(model.option_dict,SILENT_KEY,String) == SILENT_VAL

# time limit
const TIME_LIMIT = :max_wall_time
MOI.supports(::Optimizer,::MOI.TimeLimitSec) = true
MOI.set(model::Optimizer,::MOI.TimeLimitSec,value)= value isa Real ?
    MOI.set(model,MOI.RawOptimizerAttribute(String(TIME_LIMIT)),Float64(value)) : error("Invalid time limit: $value")
MOI.set(model::Optimizer,::MOI.TimeLimitSec,::Nothing)=delete!(model.option_dict,TIME_LIMIT)
MOI.get(model::Optimizer,::MOI.TimeLimitSec)=get(model.option_dict,TIME_LIMIT,Float64)

# set/get options
MOI.supports(::Optimizer,::MOI.RawOptimizerAttribute) = true
MOI.set(model::Optimizer, p::MOI.RawOptimizerAttribute, value) = (model.option_dict[Symbol(p.name)] = value)
MOI.get(model::Optimizer, p::MOI.RawOptimizerAttribute) = haskey(model.option_dict, Symbol(p.name)) ?
    (return model.option_dict[Symbol(p.name)]) : error("RawOptimizerAttribute with name $(p.name) is not set.")

# solve time
MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.ips.cnt.total_time

function MOI.empty!(model::Optimizer)
    # empty!(model.option_dict)
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
    model.nlp_dual_start = nothing
    model.ips = nothing
    model.nlp = nothing
    model.result = nothing
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) &&
        model.nlp_data.evaluator isa EmptyNLPEvaluator &&
        model.sense == MOI.FEASIBILITY_SENSE &&
        isempty(model.linear_le_constraints) &&
        isempty(model.linear_ge_constraints) &&
        isempty(model.linear_eq_constraints) &&
        isempty(model.quadratic_le_constraints) &&
        isempty(model.quadratic_ge_constraints) &&
        isempty(model.quadratic_eq_constraints)
end

function MOI.add_variable(model::Optimizer)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end
MOI.add_variables(model::Optimizer, n::Int) = [MOI.add_variable(model) for i in 1:n]

# checking inbounds
function check_inbounds(model::Optimizer, vi::MOI.VariableIndex)
    num_variables = length(model.variable_info)
    if !(1 <= vi.value <= num_variables)
        error("Invalid variable index $vi. ($num_variables variables in the model.)")
    end
end
function check_inbounds(model::Optimizer, aff::MOI.ScalarAffineFunction)
    for term in aff.terms
        check_inbounds(model, term.variable)
    end
end
function check_inbounds(model::Optimizer, quad::MOI.ScalarQuadraticFunction)
    for term in quad.affine_terms
        check_inbounds(model, term.variable)
    end
    for term in quad.quadratic_terms
        check_inbounds(model, term.variable_1)
        check_inbounds(model, term.variable_2)
    end
end

has_upper_bound(model::Optimizer, vi::MOI.VariableIndex) = model.variable_info[vi.value].has_upper_bound
has_lower_bound(model::Optimizer, vi::MOI.VariableIndex) = model.variable_info[vi.value].has_lower_bound
is_fixed(model::Optimizer, vi::MOI.VariableIndex) = model.variable_info[vi.value].is_fixed

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}})
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].upper_bound = Inf
    model.variable_info[ci.value].has_upper_bound = false
    return
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}})
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = -Inf
    model.variable_info[ci.value].has_lower_bound = false
    return
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}})
    MOI.throw_if_not_valid(model, ci)
    model.variable_info[ci.value].lower_bound = -Inf
    model.variable_info[ci.value].upper_bound = Inf
    model.variable_info[ci.value].is_fixed = false
    return
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, lt::MOI.LessThan{Float64})
    check_inbounds(model, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Float64}, typeof(lt)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.EqualTo{Float64}, typeof(lt)}(vi))
    end
    model.variable_info[vi.value].upper_bound = lt.upper
    model.variable_info[vi.value].has_upper_bound = true
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, gt::MOI.GreaterThan{Float64})
    check_inbounds(model, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Float64}, typeof(gt)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{Float64}, typeof(gt)}(vi))
    end
    model.variable_info[vi.value].lower_bound = gt.lower
    model.variable_info[vi.value].has_lower_bound = true
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, eq::MOI.EqualTo{Float64})
    check_inbounds(model, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.GreaterThan{Float64}, typeof(eq)}(vi))
    end
    if has_upper_bound(model, vi)
        throw(MOI.UpperBoundAlreadySet{MOI.LessThan{Float64}, typeof(eq)}(vi))
    end
    if is_fixed(model, vi)
        throw(MOI.LowerBoundAlreadySet{MOI.EqualTo{Float64}, typeof(eq)}(vi))
    end
    model.variable_info[vi.value].lower_bound = eq.value
    model.variable_info[vi.value].upper_bound = eq.value
    model.variable_info[vi.value].is_fixed = true
    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}}(vi.value)
end

macro define_add_constraint(function_type, set_type, prefix)
    array_name = Symbol(string(prefix) * "_constraints")
    quote
        function MOI.add_constraint(model::Optimizer, func::$function_type, set::$set_type)
            check_inbounds(model, func)
            push!(model.$(array_name), ConstraintInfo(func, set))
            return MOI.ConstraintIndex{$function_type, $set_type}(length(model.$(array_name)))
        end
    end
end

@define_add_constraint(MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64},linear_le)
@define_add_constraint(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}, linear_ge)
@define_add_constraint(MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64},linear_eq)
@define_add_constraint(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64}, quadratic_le)
@define_add_constraint(MOI.ScalarQuadraticFunction{Float64},MOI.GreaterThan{Float64}, quadratic_ge)
@define_add_constraint(MOI.ScalarQuadraticFunction{Float64},MOI.EqualTo{Float64}, quadratic_eq)

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::MOI.VariableIndex, value::Union{Real, Nothing})
    check_inbounds(model, vi)
    model.variable_info[vi.value].start = value
    return
end

function MOI.supports(model::Optimizer, ::MOI.ConstraintDualStart,
                      ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}},
                      value::Union{Real, Nothing})
    return true
end
function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}},
                 value::Union{Real, Nothing})
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    model.variable_info[vi.value].lower_bound_dual_start = value
    return
end
function MOI.supports(model::Optimizer, ::MOI.ConstraintDualStart,
                      ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}},
                      value::Union{Real, Nothing})
    return true
end
function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}},
                 value::Union{Real, Nothing})
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    model.variable_info[vi.value].upper_bound_dual_start = value
    return
end
function MOI.supports(model::Optimizer, ::MOI.ConstraintDualStart,
                      ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}},
                      value::Union{Real, Nothing})
    return true
end
function MOI.set(model::Optimizer, ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{Float64}},
                 value::Union{Real, Nothing})
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    model.variable_info[vi.value].upper_bound_dual_start = value
    model.variable_info[vi.value].lower_bound_dual_start = value
    return
end

# NLPBlock Dual start
MOI.supports(::Optimizer, ::MOI.NLPBlockDualStart) = true
MOI.set(model::Optimizer, ::MOI.NLPBlockDualStart, values) = (model.nlp_dual_start = -values)
MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData) = (model.nlp_data = nlp_data)

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction,
                 func::Union{MOI.VariableIndex, MOI.ScalarAffineFunction,MOI.ScalarQuadraticFunction})
    check_inbounds(model, func)
    model.objective = func
    return
end

linear_le_offset(model::Optimizer) = 0
linear_ge_offset(model::Optimizer) = length(model.linear_le_constraints)
linear_eq_offset(model::Optimizer) = linear_ge_offset(model) + length(model.linear_ge_constraints)
quadratic_le_offset(model::Optimizer) = linear_eq_offset(model) + length(model.linear_eq_constraints)
quadratic_ge_offset(model::Optimizer) = quadratic_le_offset(model) + length(model.quadratic_le_constraints)
quadratic_eq_offset(model::Optimizer) = quadratic_ge_offset(model) + length(model.quadratic_ge_constraints)
nlp_constraint_offset(model::Optimizer) = quadratic_eq_offset(model) + length(model.quadratic_eq_constraints)

# Convenience functions used only in optimize!
function append_to_jacobian_sparsity!(I,J,
                                      aff::MOI.ScalarAffineFunction, row, offset)
    cnt = 0
    for term in aff.terms
        I[offset+cnt]= row
        J[offset+cnt]= term.variable.value
        cnt += 1
    end
    return cnt
end

function append_to_jacobian_sparsity!(I,J,
                                      quad::MOI.ScalarQuadraticFunction, row, offset)
    cnt = 0
    for term in quad.affine_terms
        I[offset+cnt]= row
        J[offset+cnt]= term.variable.value
        cnt += 1
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        if row_idx == col_idx
            I[offset+cnt]= row
            J[offset+cnt]= row_idx.value
            cnt += 1
        else
            I[offset+cnt]= row
            J[offset+cnt]= row_idx.value
            I[offset+cnt+1]= row
            J[offset+cnt+1]= col_idx.value
            cnt += 2
        end
    end
    return cnt
end

# Refers to local variables in jacobian_structure() below.
macro append_to_jacobian_sparsity(array_name)
    escrow = esc(:row)
    escoffset = esc(:offset)
    quote
        for info in $(esc(array_name))
            $escoffset += append_to_jacobian_sparsity!($(esc(:I)), $(esc(:J)),info.func, $escrow, $escoffset)
            $escrow += 1
        end
    end
end

function jacobian_structure(model::Optimizer,I,J)
    offset = 1
    row = 1

    @append_to_jacobian_sparsity model.linear_le_constraints
    @append_to_jacobian_sparsity model.linear_ge_constraints
    @append_to_jacobian_sparsity model.linear_eq_constraints
    @append_to_jacobian_sparsity model.quadratic_le_constraints
    @append_to_jacobian_sparsity model.quadratic_ge_constraints
    @append_to_jacobian_sparsity model.quadratic_eq_constraints

    cnt = 0
    for (nlp_row, nlp_col) in MOI.jacobian_structure(model.nlp_data.evaluator)
        I[offset+cnt] = nlp_row + row - 1
        J[offset+cnt] = nlp_col
        cnt+=1
    end

    return I,J
end

append_to_hessian_sparsity!(I,J,::Union{MOI.VariableIndex,MOI.ScalarAffineFunction},offset) = 0

function append_to_hessian_sparsity!(I,J,quad::MOI.ScalarQuadraticFunction,offset)
    cnt = 0
    for term in quad.quadratic_terms
        I[offset+cnt]=term.variable_1.value
        J[offset+cnt]=term.variable_2.value
        cnt+=1
    end
    return cnt
end
function hessian_lagrangian_structure(model::Optimizer,I,J)
    offset = 1
    if !model.nlp_data.has_objective && model.objective !== nothing
        offset+=append_to_hessian_sparsity!(I,J,model.objective,offset)
    end
    for info in model.quadratic_le_constraints
        offset+=append_to_hessian_sparsity!(I,J,info.func,offset)
    end
    for info in model.quadratic_ge_constraints
        offset+=append_to_hessian_sparsity!(I,J,info.func,offset)
    end
    for info in model.quadratic_eq_constraints
        offset+=append_to_hessian_sparsity!(I,J,info.func,offset)
    end
    cnt = 0
    for (row,col) in MOI.hessian_lagrangian_structure(model.nlp_data.evaluator)
        I[offset+cnt]=row
        J[offset+cnt]=col
        cnt+=1
    end
    return I,J
end

get_nnz_hess(model::Optimizer) = get_nnz_hess_nonlinear(model) + get_nnz_hess_quadratic(model)
get_nnz_hess_nonlinear(model::Optimizer) =
    ((!model.nlp_data.has_objective && model.objective isa MOI.ScalarQuadraticFunction) ? length(model.objective.quadratic_terms) : 0) + length(MOI.hessian_lagrangian_structure(model.nlp_data.evaluator))
get_nnz_hess_quadratic(model::Optimizer) =
    (isempty(model.quadratic_eq_constraints) ? 0 : sum(length(term.func.quadratic_terms) for term in model.quadratic_eq_constraints)) +
    (isempty(model.quadratic_le_constraints) ? 0 : sum(length(term.func.quadratic_terms) for term in model.quadratic_le_constraints)) +
    (isempty(model.quadratic_ge_constraints) ? 0 : sum(length(term.func.quadratic_terms) for term in model.quadratic_ge_constraints))

get_nnz_jac(model::Optimizer) = get_nnz_jac_linear(model) + get_nnz_jac_quadratic(model) + get_nnz_jac_nonlinear(model)
get_nnz_jac_linear(model::Optimizer) =
    (isempty(model.linear_eq_constraints) ? 0 : sum(length(info.func.terms) for info in model.linear_eq_constraints)) +
    (isempty(model.linear_le_constraints) ? 0 : sum(length(info.func.terms) for info in model.linear_le_constraints)) +
    (isempty(model.linear_ge_constraints) ? 0 : sum(length(info.func.terms) for info in model.linear_ge_constraints))
get_nnz_jac_quadratic(model::Optimizer) =
    (isempty(model.quadratic_eq_constraints) ? 0 : sum(length(info.func.affine_terms) for info in model.quadratic_eq_constraints)) +
    (isempty(model.quadratic_eq_constraints) ? 0 : sum(isempty(info.func.quadratic_terms) ? 0 : sum(term.variable_1 == term.variable_2 ? 1 : 2 for term in info.func.quadratic_terms) for info in model.quadratic_eq_constraints)) +
    (isempty(model.quadratic_le_constraints) ? 0 : sum(length(info.func.affine_terms) for info in model.quadratic_le_constraints)) +
    (isempty(model.quadratic_le_constraints) ? 0 : sum(isempty(info.func.quadratic_terms) ? 0 : sum(term.variable_1 == term.variable_2 ? 1 : 2 for term in info.func.quadratic_terms) for info in model.quadratic_le_constraints)) +
    (isempty(model.quadratic_ge_constraints) ? 0 : sum(length(info.func.affine_terms) for info in model.quadratic_ge_constraints)) +
    (isempty(model.quadratic_ge_constraints) ? 0 : sum(isempty(info.func.quadratic_terms) ? 0 : sum(term.variable_1 == term.variable_2 ? 1 : 2 for term in info.func.quadratic_terms) for info in model.quadratic_ge_constraints))
get_nnz_jac_nonlinear(model::Optimizer) =
    (isempty(model.nlp_data.constraint_bounds) ? 0 : length(MOI.jacobian_structure(model.nlp_data.evaluator)))

function eval_function(var::MOI.VariableIndex, x)
    return x[var.value]
end

function eval_function(aff::MOI.ScalarAffineFunction, x)
    function_value = aff.constant
    for term in aff.terms
        function_value += term.coefficient*x[term.variable.value]
    end
    return function_value
end

function eval_function(quad::MOI.ScalarQuadraticFunction, x)
    function_value = quad.constant

    for term in quad.affine_terms
        function_value += term.coefficient*x[term.variable.value]
    end

    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        coefficient = term.coefficient
        if row_idx == col_idx
            function_value += 0.5*coefficient*x[row_idx.value]*x[col_idx.value]
        else
            function_value += coefficient*x[row_idx.value]*x[col_idx.value]
        end
    end
    return function_value
end

function eval_objective(model::Optimizer, x)
    # The order of the conditions is important. NLP objectives override regular
    # objectives.
    if model.nlp_data.has_objective
        return MOI.eval_objective(model.nlp_data.evaluator, x)
    elseif model.objective !== nothing
        return eval_function(model.objective, x)
    else
        # No objective function set. This could happen with FEASIBILITY_SENSE.
        return 0.0
    end
end

function fill_gradient!(grad, x, var::MOI.VariableIndex)
    fill!(grad, 0.0)
    grad[var.value] = 1.0
end

function fill_gradient!(grad, x, aff::MOI.ScalarAffineFunction{Float64})
    fill!(grad, 0.0)
    for term in aff.terms
        grad[term.variable.value] += term.coefficient
    end
end

function fill_gradient!(grad, x, quad::MOI.ScalarQuadraticFunction{Float64})
    fill!(grad, 0.0)
    for term in quad.affine_terms
        grad[term.variable.value] += term.coefficient
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        coefficient = term.coefficient
        if row_idx == col_idx
            grad[row_idx.value] += coefficient*x[row_idx.value]
        else
            grad[row_idx.value] += coefficient*x[col_idx.value]
            grad[col_idx.value] += coefficient*x[row_idx.value]
        end
    end
end

function eval_objective_gradient(model::Optimizer, grad, x)
    if model.nlp_data.has_objective
        MOI.eval_objective_gradient(model.nlp_data.evaluator, grad, x)
    elseif model.objective !== nothing
        fill_gradient!(grad, x, model.objective)
    else
        fill!(grad, 0.0)
    end
    return
end

# Refers to local variables in eval_constraint() below.
macro eval_function(array_name)
    escrow = esc(:row)
    quote
        for info in $(esc(array_name))
            $(esc(:g))[$escrow] = eval_function(info.func, $(esc(:x)))
            $escrow += 1
        end
    end
end
function eval_constraint(model::Optimizer, g, x)
    row = 1
    @eval_function model.linear_le_constraints
    @eval_function model.linear_ge_constraints
    @eval_function model.linear_eq_constraints
    @eval_function model.quadratic_le_constraints
    @eval_function model.quadratic_ge_constraints
    @eval_function model.quadratic_eq_constraints
    nlp_g = view(g, row:length(g))
    MOI.eval_constraint(model.nlp_data.evaluator, nlp_g, x)
    return
end

function fill_constraint_jacobian!(values, start_offset, x, aff::MOI.ScalarAffineFunction)
    num_coefficients = length(aff.terms)
    for i in 1:num_coefficients
        values[start_offset+i] = aff.terms[i].coefficient
    end
    return num_coefficients
end

function fill_constraint_jacobian!(values, start_offset, x, quad::MOI.ScalarQuadraticFunction)
    num_affine_coefficients = length(quad.affine_terms)
    for i in 1:num_affine_coefficients
        values[start_offset+i] = quad.affine_terms[i].coefficient
    end
    num_quadratic_coefficients = 0
    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        coefficient = term.coefficient
        if row_idx == col_idx
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+1] = coefficient*x[col_idx.value]
            num_quadratic_coefficients += 1
        else
            # Note that the order matches the Jacobian sparsity pattern.
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+1] = coefficient*x[col_idx.value]
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+2] = coefficient*x[row_idx.value]
            num_quadratic_coefficients += 2
        end
    end
    return num_affine_coefficients + num_quadratic_coefficients
end

# Refers to local variables in eval_constraint_jacobian() below.
macro fill_constraint_jacobian(array_name)
    esc_offset = esc(:offset)
    quote
        for info in $(esc(array_name))
            $esc_offset += fill_constraint_jacobian!($(esc(:values)),
                                                     $esc_offset, $(esc(:x)),
                                                     info.func)
        end
    end
end

function eval_constraint_jacobian(model::Optimizer, values, x)
    offset = 0
    @fill_constraint_jacobian model.linear_le_constraints
    @fill_constraint_jacobian model.linear_ge_constraints
    @fill_constraint_jacobian model.linear_eq_constraints
    @fill_constraint_jacobian model.quadratic_le_constraints
    @fill_constraint_jacobian model.quadratic_ge_constraints
    @fill_constraint_jacobian model.quadratic_eq_constraints

    nlp_values = view(values, 1+offset:length(values))
    MOI.eval_constraint_jacobian(model.nlp_data.evaluator, nlp_values, x)
    return
end

function fill_hessian_lagrangian!(values, start_offset, scale_factor,
                                  ::Union{MOI.VariableIndex,
                                          MOI.ScalarAffineFunction,Nothing})
    return 0
end

function fill_hessian_lagrangian!(values, start_offset, scale_factor,
                                  quad::MOI.ScalarQuadraticFunction)
    for i in 1:length(quad.quadratic_terms)
        values[start_offset + i] = scale_factor*quad.quadratic_terms[i].coefficient
    end
    return length(quad.quadratic_terms)
end

function eval_hessian_lagrangian(model::Optimizer, values, x, obj_factor, lambda)
    offset = 0
    if !model.nlp_data.has_objective
        offset += fill_hessian_lagrangian!(values, 0, obj_factor,
                                           model.objective)
    end
    for (i, info) in enumerate(model.quadratic_le_constraints)
        offset += fill_hessian_lagrangian!(values, offset, lambda[i+quadratic_le_offset(model)], info.func)
    end
    for (i, info) in enumerate(model.quadratic_ge_constraints)
        offset += fill_hessian_lagrangian!(values, offset, lambda[i+quadratic_ge_offset(model)], info.func)
    end
    for (i, info) in enumerate(model.quadratic_eq_constraints)
        offset += fill_hessian_lagrangian!(values, offset, lambda[i+quadratic_eq_offset(model)], info.func)
    end
    nlp_values = view(values, 1 + offset : length(values))
    nlp_lambda = view(lambda, 1 + nlp_constraint_offset(model) : length(lambda))
    MOI.eval_hessian_lagrangian(model.nlp_data.evaluator, nlp_values, x, obj_factor, nlp_lambda)
end

MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.result === nothing ?
    MOI.OPTIMIZE_NOT_CALLED : termination_status(model.result)
MOI.get(model::Optimizer, ::MOI.RawStatusString) = string(model.result.status)
MOI.get(model::Optimizer, ::MOI.ResultCount) = (model.result !== nothing) ? 1 : 0
MOI.get(model::Optimizer, attr::MOI.PrimalStatus) = !(1 <= attr.result_index <= MOI.get(model, MOI.ResultCount())) ?
    MOI.NO_SOLUTION : primal_status(model.result)
MOI.get(model::Optimizer, attr::MOI.DualStatus) = !(1 <= attr.result_index <= MOI.get(model, MOI.ResultCount())) ?
    MOI.NO_SOLUTION : dual_status(model.result)

const status_moi_dict = Dict(
    SOLVE_SUCCEEDED => MOI.LOCALLY_SOLVED,
    SOLVED_TO_ACCEPTABLE_LEVEL => MOI.ALMOST_LOCALLY_SOLVED,
    SEARCH_DIRECTION_BECOMES_TOO_SMALL => MOI.SLOW_PROGRESS,
    DIVERGING_ITERATES => MOI.INFEASIBLE_OR_UNBOUNDED,
    INFEASIBLE_PROBLEM_DETECTED => MOI.LOCALLY_INFEASIBLE,
    MAXIMUM_ITERATIONS_EXCEEDED => MOI.ITERATION_LIMIT,
    MAXIMUM_WALLTIME_EXCEEDED => MOI.TIME_LIMIT,
    INITIAL => MOI.OPTIMIZE_NOT_CALLED,
    RESTORATION_FAILED => MOI.NUMERICAL_ERROR,
    INVALID_NUMBER_DETECTED => MOI.INVALID_MODEL,
    ERROR_IN_STEP_COMPUTATION => MOI.NUMERICAL_ERROR,
    NOT_ENOUGH_DEGREES_OF_FREEDOM => MOI.INVALID_MODEL,
    USER_REQUESTED_STOP => MOI.INTERRUPTED,
    INTERNAL_ERROR => MOI.OTHER_ERROR)
const status_primal_dict = Dict(
    SOLVE_SUCCEEDED => MOI.FEASIBLE_POINT,
    SOLVED_TO_ACCEPTABLE_LEVEL => MOI.NEARLY_FEASIBLE_POINT,
    INFEASIBLE_PROBLEM_DETECTED => MOI.INFEASIBLE_POINT)
const status_dual_dict = Dict(
    SOLVE_SUCCEEDED => MOI.FEASIBLE_POINT,
    SOLVED_TO_ACCEPTABLE_LEVEL => MOI.NEARLY_FEASIBLE_POINT,
    INFEASIBLE_PROBLEM_DETECTED => MOI.INFEASIBLE_POINT)

termination_status(result::MadNLPExecutionStats) = haskey(status_moi_dict,result.status) ? status_moi_dict[result.status] : MOI.UNKNOWN_RESULT_STATUS
termination_status(ips::InteriorPointSolver) = haskey(status_moi_dict,ips.status) ? status_moi_dict[ips.status] : MOI.UNKNOWN_RESULT_STATUS
primal_status(result::MadNLPExecutionStats) = haskey(status_primal_dict,result.status) ?
    status_primal_dict[result.status] : MOI.UNKNOWN_RESULT_STATUS
dual_status(result::MadNLPExecutionStats) = haskey(status_dual_dict,result.status) ?
    status_dual_dict[result.status] : MOI.UNKNOWN_RESULT_STATUS

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    scale = (model.sense == MOI.MAX_SENSE) ? -1 : 1
    return scale * model.result.objective
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(model, attr)
    check_inbounds(model, vi)
    return model.result.solution[vi.value]
end

macro define_constraint_primal(function_type, set_type, prefix)
    constraint_array = Symbol(string(prefix) * "_constraints")
    offset_function = Symbol(string(prefix) * "_offset")
    quote
        function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                         ci::MOI.ConstraintIndex{$function_type, $set_type})
            MOI.check_result_index_bounds(model, attr)
            if !(1 <= ci.value <= length(model.$(constraint_array)))
                error("Invalid constraint index ", ci.value)
            end
            return model.result.constraints[ci.value + $offset_function(model)]
        end
    end
end

@define_constraint_primal(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}, linear_le)
@define_constraint_primal(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}, linear_ge)
@define_constraint_primal(MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}, linear_eq)
@define_constraint_primal(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64}, quadratic_le)
@define_constraint_primal(MOI.ScalarQuadraticFunction{Float64},MOI.GreaterThan{Float64}, quadratic_ge)
@define_constraint_primal(MOI.ScalarQuadraticFunction{Float64},MOI.EqualTo{Float64}, quadratic_eq)

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,
                                         MOI.LessThan{Float64}})
    MOI.check_result_index_bounds(model, attr)
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !has_upper_bound(model, vi)
        error("Variable $vi has no upper bound -- ConstraintPrimal not defined.")
    end
    return model.result.solution[vi.value]
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,
                                         MOI.GreaterThan{Float64}})
    MOI.check_result_index_bounds(model, attr)
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !has_lower_bound(model, vi)
        error("Variable $vi has no lower bound -- ConstraintPrimal not defined.")
    end
    return model.result.solution[vi.value]
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,
                                         MOI.EqualTo{Float64}})
    MOI.check_result_index_bounds(model, attr)
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !is_fixed(model, vi)
        error("Variable $vi is not fixed -- ConstraintPrimal not defined.")
    end
    return model.result.solution[vi.value]
end

macro define_constraint_dual(function_type, set_type, prefix)
    constraint_array = Symbol(string(prefix) * "_constraints")
    offset_function = Symbol(string(prefix) * "_offset")
    quote
        function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
                         ci::MOI.ConstraintIndex{$function_type, $set_type})
            MOI.check_result_index_bounds(model, attr)
            if !(1 <= ci.value <= length(model.$(constraint_array)))
                error("Invalid constraint index ", ci.value)
            end
            return -1 * model.result.multipliers[ci.value + $offset_function(model)]
        end
    end
end
@define_constraint_dual(MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}, linear_le)
@define_constraint_dual(MOI.ScalarAffineFunction{Float64},MOI.GreaterThan{Float64}, linear_ge)
@define_constraint_dual(MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}, linear_eq)
@define_constraint_dual(MOI.ScalarQuadraticFunction{Float64},MOI.LessThan{Float64}, quadratic_le)
@define_constraint_dual(MOI.ScalarQuadraticFunction{Float64},MOI.GreaterThan{Float64}, quadratic_ge)
@define_constraint_dual(MOI.ScalarQuadraticFunction{Float64},MOI.EqualTo{Float64}, quadratic_eq)

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}})
    MOI.check_result_index_bounds(model, attr)
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    has_upper_bound(model, vi) || error("Variable $vi has no upper bound -- ConstraintDual not defined.")
    return -1 * model.result.multipliers_U[vi.value] # MOI convention is for feasible LessThan duals to be nonpositive.
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,
                                         MOI.GreaterThan{Float64}})
    MOI.check_result_index_bounds(model, attr)
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    has_lower_bound(model, vi) || error("Variable $vi has no lower bound -- ConstraintDual not defined.")
    return model.result.multipliers_L[vi.value]
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,
                                         MOI.EqualTo{Float64}})
    MOI.check_result_index_bounds(model, attr)
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    if !is_fixed(model, vi)
        error("Variable $vi is not fixed -- ConstraintDual not defined.")
    end
    return model.result.multipliers_L[vi.value] - model.result.multipliers_U[vi.value]
end

function MOI.get(model::Optimizer, attr::MOI.NLPBlockDual)
    MOI.check_result_index_bounds(model, attr)
    return -1 * model.result.multipliers[(1 + nlp_constraint_offset(model)):end]
end

function num_constraints(model::Optimizer)
    m_linear_le = length(model.linear_le_constraints)
    m_linear_ge = length(model.linear_ge_constraints)
    m_linear_eq = length(model.linear_eq_constraints)
    m_quadratic_le = length(model.quadratic_le_constraints)
    m_quadratic_ge = length(model.quadratic_ge_constraints)
    m_quadratic_eq = length(model.quadratic_eq_constraints)
    m_nl = length(model.nlp_data.constraint_bounds)
    return m_linear_le + m_linear_ge + m_linear_eq + m_quadratic_le + m_quadratic_ge + m_quadratic_eq + m_nl
end

function is_jac_hess_constant(model::Optimizer)
    isempty(model.nlp_data.constraint_bounds) || return (false,false)
    isempty(model.quadratic_le_constraints) || return (false,false)
    isempty(model.quadratic_ge_constraints) || return (false,false)
    isempty(model.quadratic_eq_constraints) || return (false,false)
    return (true,model.nlp_data.has_objective ? false : true)
end

zero_if_nothing(x) = x == nothing ? 0. : x

function set_x!(model,x,xl,xu)

    for i=1:length(model.variable_info)
        info = model.variable_info[i]
        x[i]  = info.start == nothing ? 0. : info.start
        xl[i] = info.lower_bound
        xu[i] = info.upper_bound
        i += 1
    end
end

function set_g!(model::Optimizer,l,gl,gu)
    i = 1
    for info in model.linear_le_constraints
        l[i] = zero_if_nothing(info.dual_start)
        gl[i]= -Inf
        gu[i]= info.set.upper
        i+=1
    end
    for info in model.linear_ge_constraints
        l[i] = zero_if_nothing(info.dual_start)
        gl[i]= info.set.lower
        gu[i]= Inf
        i+=1
    end
    for info in model.linear_eq_constraints
        l[i] = zero_if_nothing(info.dual_start)
        gl[i]= info.set.value
        gu[i]= info.set.value
        i+=1
    end
    for info in model.quadratic_le_constraints
        l[i] = zero_if_nothing(info.dual_start)
        gl[i]= -Inf
        gu[i]= info.set.upper
        i+=1
    end
    for info in model.quadratic_ge_constraints
        l[i] = zero_if_nothing(info.dual_start)
        gl[i]= info.set.lower
        gu[i]= Inf
        i+=1
    end
    for info in model.quadratic_eq_constraints
        l[i] = zero_if_nothing(info.dual_start)
        gl[i]= info.set.value
        gu[i]= info.set.value
        i+=1
    end
    j=0
    for bound in model.nlp_data.constraint_bounds
        gl[i+j]= bound.lower
        gu[i+j]= bound.upper
        j+=1
    end
    k=0
    if model.nlp_dual_start != nothing
        for v in model.nlp_dual_start
            l[i+k] = v
            k+=1
        end
    else
        l[i:end] .= 0.
    end
end

num_variables(model::Optimizer) = length(model.variable_info)

struct MOIModel{T} <: AbstractNLPModel{T,Vector{T}}
    meta::NLPModelMeta{T, Vector{T}}
    model::Optimizer
    counters::NLPModelsCounters
end

obj(nlp::MOIModel,x::StrideOneVector{Float64}) = eval_objective(nlp.model,x)
grad!(nlp::MOIModel,x::StrideOneVector{Float64},f::StrideOneVector{Float64}) =
    eval_objective_gradient(nlp.model,f,x)
cons!(nlp::MOIModel,x::StrideOneVector{Float64},c::StrideOneVector{Float64}) =
    eval_constraint(nlp.model,c,x)
jac_coord!(nlp::MOIModel,x::StrideOneVector{Float64},jac::StrideOneVector{Float64})=
    eval_constraint_jacobian(nlp.model,jac,x)
hess_coord!(nlp::MOIModel,x::StrideOneVector{Float64},l::StrideOneVector{Float64},hess::StrideOneVector{Float64};
            obj_weight::Float64=1.) = eval_hessian_lagrangian(nlp.model,hess,x,obj_weight,l)
function hess_structure!(nlp::MOIModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
    return hessian_lagrangian_structure(nlp.model,I,J)
end
function jac_structure!(nlp::MOIModel, I::AbstractVector{T}, J::AbstractVector{T}) where T
    return jacobian_structure(nlp.model,I,J)
end

function MOIModel(model::Optimizer)
    :Hess in MOI.features_available(model.nlp_data.evaluator) || error("Hessian information is needed.")
    MOI.initialize(model.nlp_data.evaluator, [:Grad,:Hess,:Jac])

    nvar = num_variables(model)
    ncon = num_constraints(model)
    nnzh = get_nnz_hess(model)
    nnzj = get_nnz_jac(model)

    x0  = Vector{Float64}(undef,nvar)
    xl = Vector{Float64}(undef,nvar)
    xu = Vector{Float64}(undef,nvar)

    y0 = Vector{Float64}(undef,ncon)
    gl = Vector{Float64}(undef,ncon)
    gu = Vector{Float64}(undef,ncon)

    set_x!(model,x0,xl,xu)
    set_g!(model,y0,gl,gu)

    model.option_dict[:jacobian_constant], model.option_dict[:hessian_constant] = is_jac_hess_constant(model)
    model.option_dict[:dual_initialized] = !iszero(y0)

    return MOIModel(
        NLPModelMeta(
            nvar,
            x0 = x0,
            lvar = xl,
            uvar = xu,
            ncon = ncon,
            y0 = y0,
            lcon = gl,
            ucon = gu,
            nnzj = nnzj,
            nnzh = nnzh,
            minimize = model.sense == MOI.MIN_SENSE
        ),
        model,NLPModelsCounters())
end

function MOI.optimize!(model::Optimizer)

    model.nlp = MOIModel(model)
    model.ips = InteriorPointSolver(model.nlp;option_dict=copy(model.option_dict))
    model.result = optimize!(model.ips)
    model.solve_time = model.ips.cnt.total_time

    return
end


text:
# ---
# title: 955. Delete Columns to Make Sorted II
# id: problem955
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/delete-columns-to-make-sorted-ii/description/>
# hidden: true
# ---
#
# We are given an array `A` of `N` lowercase letter strings, all of the same
# length.
#
# Now, we may choose any set of deletion indices, and for each string, we delete
# all the characters in those indices.
#
# For example, if we have an array `A = ["abcdef","uvwxyz"]` and deletion
# indices `{0, 2, 3}`, then the final array after deletions is `["bef","vyz"]`.
#
# Suppose we chose a set of deletion indices `D` such that after deletions, the
# final array has its elements in **lexicographic** order (`A[0] <= A[1] <= A[2]
# ... <= A[A.length - 1]`).
#
# Return the minimum possible value of `D.length`.
#
#
#
# **Example 1:**
#
#
#
#     Input: ["ca","bb","ac"]
#     Output: 1
#     Explanation:
#     After deleting the first column, A = ["a", "b", "c"].
#     Now A is in lexicographic order (ie. A[0] <= A[1] <= A[2]).
#     We require at least 1 deletion since initially A was not in lexicographic order, so the answer is 1.
#
#
# **Example 2:**
#
#
#
#     Input: ["xc","yb","za"]
#     Output: 0
#     Explanation:
#     A is already in lexicographic order, so we don't need to delete anything.
#     Note that the rows of A are not necessarily in lexicographic order:
#     ie. it is NOT necessarily true that (A[0][0] <= A[0][1] <= ...)
#
#
# **Example 3:**
#
#
#
#     Input: ["zyx","wvu","tsr"]
#     Output: 3
#     Explanation:
#     We have to delete every column.
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 100`
#   2. `1 <= A[i].length <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module TimeSpans

using Dates

export TimeSpan, start, stop, istimespan, translate, overlaps,
       shortest_timespan_containing, duration, index_from_time,
       time_from_index

const NS_IN_SEC = Dates.value(Nanosecond(Second(1)))  # Number of nanoseconds in one second

#####
##### `TimeSpan`
#####

"""
    TimeSpan(start, stop)

Return `TimeSpan(Nanosecond(start), Nanosecond(stop))` representing the interval `[start, stop)`.

If `start == stop`, a single `Nanosecond` is added to `stop` since `stop` is an exclusive
upper bound and TimeSpan operations only generally support up to nanosecond precision anyway.

The benefit of this type over e.g. `Nanosecond(start):Nanosecond(1):Nanosecond(stop)` is
that instances of this type are guaranteed to obey `TimeSpans.start(x) < TimeSpans.stop(x)`
by construction.
"""
struct TimeSpan
    start::Nanosecond
    stop::Nanosecond
    function TimeSpan(start::Nanosecond, stop::Nanosecond)
        stop += Nanosecond(start == stop)
        start < stop || throw(ArgumentError("start(span) < stop(span) must be true, got $start and $stop"))
        return new(start, stop)
    end
    TimeSpan(start, stop) = TimeSpan(Nanosecond(start), Nanosecond(stop))
end

"""
    TimeSpan(x)

Return `TimeSpan(start(x), stop(x))`.
"""
TimeSpan(x) = TimeSpan(start(x), stop(x))

Base.in(x::TimePeriod, y::TimeSpan) = start(y) <= x < stop(y)

# work around <https://github.com/JuliaLang/julia/issues/40311>:
Base.findall(pred::Base.Fix2{typeof(in), TimeSpan}, obj::Union{Tuple, AbstractArray}) = invoke(findall, Tuple{Function, typeof(obj)}, pred, obj)

#####
##### pretty printing
#####

function nanosecond_to_periods(ns::Integer)
    μs, ns = divrem(ns, 1000)
    ms, μs = divrem(μs, 1000)
    s, ms = divrem(ms, 1000)
    m, s = divrem(s, 60)
    hr, m = divrem(m, 60)
    return (hr, m, s, ms, μs, ns)
end

format_duration(t::Period) = format_duration(convert(Nanosecond, t).value)

function format_duration(ns::Integer)
    sig = signbit(ns) ? "-" : ""
    hr, m, s, ms, μs, ns = nanosecond_to_periods(abs(ns))
    hr = lpad(hr, 2, '0')
    m = lpad(m, 2, '0')
    s = lpad(s, 2, '0')
    ms = lpad(ms, 3, '0')
    μs = lpad(μs, 3, '0')
    ns = lpad(ns, 3, '0')
    return string(sig, hr, ':', m, ':', s, '.', ms, μs, ns)
end

function Base.show(io::IO, w::TimeSpan)
    start_string = format_duration(start(w))
    stop_string = format_duration(stop(w))
    return print(io, "TimeSpan(", start_string, ", ", stop_string, ')')
end

#####
##### generic TimeSpans.jl interface
#####

"""
    istimespan(x)

Return `true` if `x` has been declared to support `TimeSpans.start(x)` and `TimeSpans.stop(x)`,
return `false` otherwise.

Types that overload `TimeSpans.start`/`TimeSpans.stop` should also overload `istimespan`.
"""
istimespan(::Any) = false
istimespan(::TimeSpan) = true
istimespan(::Period) = true

"""
    start(span)

Return the inclusive lower bound of `span` as a `Nanosecond` value.
"""
start(span::TimeSpan) = span.start
start(t::Period) = convert(Nanosecond, t)

"""
    stop(span)

Return the exclusive upper bound of `span` as a `Nanosecond` value.
"""
stop(span::TimeSpan) = span.stop
stop(t::Period) = convert(Nanosecond, t) + Nanosecond(1)

#####
##### generic utilities
#####

"""
    translate(span, by::Period)

Return `TimeSpan(start(span) + by, stop(span) + by)`.
"""
function translate(span, by::Period)
    by = convert(Nanosecond, by)
    return TimeSpan(start(span) + by, stop(span) + by)
end

"""
    TimeSpans.contains(a, b)

Return `true` if the timespan `b` lies entirely within the timespan `a`, return `false` otherwise.
"""
contains(a, b) = start(a) <= start(b) && stop(a) >= stop(b)

"""
    overlaps(a, b)

Return `true` if the timespan `a` and the timespan `b` overlap, return `false` otherwise.
"""
function overlaps(a, b)
    starts_earlier, starts_later = ifelse(start(b) > start(a), (a, b), (b, a))
    return stop(starts_earlier) > start(starts_later)
end

"""
    shortest_timespan_containing(spans)

Return the shortest possible `TimeSpan` containing all timespans in `spans`.

`spans` is assumed to be an iterable of timespans.
"""
function shortest_timespan_containing(spans)
    isempty(spans) && throw(ArgumentError("input iterator must be nonempty"))
    lo, hi = Nanosecond(typemax(Int64)), Nanosecond(0)
    for span in spans
        lo = min(start(span), lo)
        hi = max(stop(span), hi)
    end
    return TimeSpan(lo, hi)
end

"""
    duration(span)

Return `stop(span) - start(span)`.
"""
duration(span) = stop(span) - start(span)

"""
    TimeSpans.nanoseconds_per_sample(sample_rate)

Given `sample_rate` in Hz, return the number of nanoseconds corresponding to one sample.
"""
nanoseconds_per_sample(sample_rate) = inv(sample_rate) * NS_IN_SEC

"""
    index_from_time(sample_rate, sample_time::Period)

Given `sample_rate` in Hz, return the integer index of the most recent sample
taken at `sample_time`. Note that `sample_time` must be non-negative and support
`convert(Nanosecond, sample_time)`.

Examples:

```jldoctest
julia> index_from_time(1, Second(0))
1

julia> index_from_time(1, Second(1))
2

julia> index_from_time(100, Millisecond(999))
100

julia> index_from_time(100, Millisecond(1000))
101
```
"""
function index_from_time(sample_rate, sample_time::Period)
    time_in_nanoseconds = convert(Nanosecond, sample_time).value
    time_in_nanoseconds >= 0 || throw(ArgumentError("`sample_time` must be >= 0 nanoseconds"))
    time_in_seconds = time_in_nanoseconds / NS_IN_SEC
    return floor(Int, time_in_seconds * sample_rate) + 1
end

"""
    index_from_time(sample_rate, span)

Return the `UnitRange` of indices corresponding to `span` given `sample_rate` in Hz:

```jldoctest
julia> index_from_time(100, TimeSpan(Second(0), Second(1)))
1:100

julia> index_from_time(100, TimeSpan(Second(1)))
101:101

julia> index_from_time(100, TimeSpan(Second(3), Second(6)))
301:600
```
"""
function index_from_time(sample_rate, span)
    i = index_from_time(sample_rate, start(span))
    j = index_from_time(sample_rate, stop(span))
    j = i == j ? j : (j - 1)
    return i:j
end

"""
    time_from_index(sample_rate, sample_index)

Given `sample_rate` in Hz and assuming `sample_index > 0`, return the earliest
`Nanosecond` containing `sample_index`.

Examples:

```jldoctest
julia> time_from_index(1, 1)
0 nanoseconds

julia> time_from_index(1, 2)
1000000000 nanoseconds

julia> time_from_index(100, 100)
990000000 nanoseconds

julia> time_from_index(100, 101)
1000000000 nanoseconds
```
"""
function time_from_index(sample_rate, sample_index)
    sample_index > 0 || throw(ArgumentError("`sample_index` must be > 0"))
    return Nanosecond(ceil(Int, (sample_index - 1) * nanoseconds_per_sample(sample_rate)))
end

"""
    time_from_index(sample_rate, sample_range::AbstractUnitRange)

Return the `TimeSpan` corresponding to `sample_range` given `sample_rate` in Hz:

```jldoctest
julia> time_from_index(100, 1:100)
TimeSpan(0 nanoseconds, 1000000000 nanoseconds)

julia> time_from_index(100, 101:101)
TimeSpan(1000000000 nanoseconds, 1000000001 nanoseconds)

julia> time_from_index(100, 301:600)
TimeSpan(3000000000 nanoseconds, 6000000000 nanoseconds)
```
"""
function time_from_index(sample_rate, sample_range::AbstractUnitRange)
    i, j = first(sample_range), last(sample_range)
    j = j == i ? j : j + 1
    return TimeSpan(time_from_index(sample_rate, i),
                    time_from_index(sample_rate, j))
end

end # module


text:
## LSMR works on any matrix with the following methods:
## eltype, size, A_mul_B!, Ac_mul_B!
## LSMR works on any vector with the following methods:
## eltype, length, scale!, norm

#############################################################################
##
## Define preconditioned matrix
##
#############################################################################

type PreconditionedMatrix{TA, Tp, Tx}
    A::TA
    preconditioner::Tp
    tmp::Tx
    tmp2::Tx
end
eltype(A::PreconditionedMatrix) = eltype(A.A)
size(A::PreconditionedMatrix, i::Integer) = size(A.A, i)
function A_mul_B!{TA, Tp, Tx}(α::Number, pm::PreconditionedMatrix{TA, Tp, Tx}, a::Tx,
                β::Number, b)
    A_ldiv_B!(a, pm.preconditioner, pm.tmp)
    A_mul_B!(α, pm.A, pm.tmp, β, b)
    return b
end
function Ac_mul_B!{TA, Tp, Tx}(α::Number, pm::PreconditionedMatrix{TA, Tp, Tx}, a,
                β::Number, b::Tx)
    T = eltype(b)
    β = convert(T, β)
    Ac_mul_B!(one(T), pm.A, a, zero(T), pm.tmp)
    A_ldiv_B!(pm.tmp, pm.preconditioner, pm.tmp2)
    if β != one(T)
        if β == zero(T)
            fill!(b, β)
        else
            scale!(b, β)
        end
    end
    axpy!(α, pm.tmp2, b)
    return b
end

#############################################################################
##
## Define dampened matrix
##
#############################################################################

type DampenedVector{Ty, Tx}
    y::Ty # dimension of f(x)
    x::Tx # dimension of x
end
eltype(a::DampenedVector) =  promote_type(eltype(a.y), eltype(a.x))
length(a::DampenedVector) = length(a.y) + length(a.x)
function scale!(a::DampenedVector, α::Number)
    scale!(a.y, α)
    scale!(a.x, α)
    return a
end
norm(a::DampenedVector) = sqrt(norm(a.y)^2 + norm(a.x)^2)

type DampenedMatrix{TA, Tx}
    A::TA
    diagonal::Tx
end
eltype(A::DampenedMatrix) = promote_type(eltype(A.A), eltype(A.diagonal))
function size(A::DampenedMatrix, dim::Integer)
    m, n = size(A.A)
    l = length(A.diagonal)
    dim == 1 ? (m + l) :
    dim == 2 ? n : 1
end
function A_mul_B!{TA, Tx, Ty}(α::Number, mw::DampenedMatrix{TA, Tx}, a::Tx,
                β::Number, b::DampenedVector{Ty, Tx})
    if β != 1
        scale!(b, β)
    end
    A_mul_B!(α, mw.A, a, 1, b.y)
    map!((z, x, y)-> z + α * x * y, b.x, b.x, a, mw.diagonal)
    return b
end
function Ac_mul_B!{TA, Tx, Ty}(α::Number, mw::DampenedMatrix{TA, Tx}, a::DampenedVector{Ty, Tx},
                β::Number, b::Tx)
    T = eltype(b)
    β = convert(T, β)
    if β != one(T)
        if β == zero(T)
            fill!(b, β)
        else
            scale!(b, β)
        end
    end
    Ac_mul_B!(α, mw.A, a.y, one(T), b)
    map!((z, x, y)-> z + α * x * y, b, b, a.x, mw.diagonal)
    return b
end

#############################################################################
##
## Define inverse preconditioner
##
##############################################################################

type InverseDiagonal{Tx}
    _::Tx
end
function A_ldiv_B!(x, ID::InverseDiagonal, y)
    map!(*, y, x, ID._)
end

#############################################################################
##
## solve J'J \ J'y
##
## we use LSMR on Ax = y
## with A = J / sqrt(diag(J'J)) (diagonal preconditioner)
##
## LSMR works on any matrix with the following methods:
## eltype, size, A_mul_B!, Ac_mul_B!
##
##############################################################################
function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg, solver::LSMR{Void, Void})
    preconditioner! = function(x, J, out)
        colsumabs2!(out._, J)
        Tout = eltype(out._)
        map!(x -> x > zero(Tout) ? 1 / sqrt(x) : zero(Tout), out._, out._)
        return out
    end
    preconditioner = InverseDiagonal(_zeros(nls.x))
    return preconditioner!, preconditioner
end
function getpreconditioner(nls::LeastSquaresProblem, optimizer::Dogleg, solver::LSMR)
    return solver.preconditioner!, solver.preconditioner
end

type LSMRAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Ty} <: AbstractAllocatedSolver
    preconditioner!::Tx0
    preconditioner::Tx1
    tmp::Tx2
    tmp2::Tx22
    v::Tx3
    h::Tx4
    hbar::Tx5
    u::Ty
    function LSMRAllocatedSolver(preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, u)
        new(preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, u)
    end
end

function LSMRAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Ty}(preconditioner!::Tx0, preconditioner::Tx1, tmp::Tx2, tmp2::Tx22, v::Tx3, h::Tx4, hbar::Tx5, u::Ty)
    LSMRAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Ty}(preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, u)
end

function AbstractAllocatedSolver(nls::LeastSquaresProblem, optimizer::Dogleg, solver::LSMR)
    preconditioner!, preconditioner = getpreconditioner(nls, optimizer, solver)
    LSMRAllocatedSolver(preconditioner!, preconditioner, _zeros(nls.x), _zeros(nls.x),
        _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.y))
end

function A_ldiv_B!(x, J, y, A::LSMRAllocatedSolver)
    preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, u = A.preconditioner!, A.preconditioner, A.tmp, A.tmp2, A.v, A.h, A.hbar, A.u

    # prepare x
    fill!(x, 0)

    # prepare b
    copy!(u, y)

    # prepare A
    fill!(tmp, 0)
    preconditioner!(x, J, preconditioner)
    A = PreconditionedMatrix(J, preconditioner, tmp, tmp2)

    # solve
    x, ch = lsmr!(x, A, u, v, h, hbar)
    A_ldiv_B!(x, preconditioner, tmp)
    copy!(x, tmp)
    return x, ch.mvps
end

##############################################################################
##
## solve (J'J + diagm(damp)) \ J'y (used in LevenbergMarquardt)
## No need to solve exactly :
## "An Inexact Levenberg-Marquardt Method for Large Sparse Nonlinear Least Squares"
## Weight Holt (1985)
##
## We use LSMR on A x = b with
## A = |J          |  + diagonal preconditioner
##     |diag(damp) |
## b = vcat(y, zeros(damp))
##

##
##############################################################################
function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt, ::LSMR{Void, Void})
    preconditioner! = function(x, J, damp, out)
        colsumabs2!(out._, J)
        Tout = eltype(out._)
        axpy!(one(Tout), damp, out._)
        map!(x -> x > zero(Tout) ? 1 / sqrt(x) : zero(Tout), out._, out._)
        return out
    end
    preconditioner = InverseDiagonal(_zeros(nls.x))
    return preconditioner!, preconditioner
end

function getpreconditioner(nls::LeastSquaresProblem, optimizer::LevenbergMarquardt, solver::LSMR)
    return solver.preconditioner!, solver.preconditioner
end

type LSMRDampenedAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Tx6, Ty} <: AbstractAllocatedSolver
    preconditioner!::Tx0
    preconditioner::Tx1
    tmp::Tx2
    tmp2::Tx22
    v::Tx3
    h::Tx4
    hbar::Tx5
    zerosvector::Tx6
    u::Ty
    function LSMRDampenedAllocatedSolver(preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, zerosvector, u)
        new(preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, zerosvector, u)
    end
end
function LSMRDampenedAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Tx6, Ty}(preconditioner!::Tx0, preconditioner::Tx1, tmp::Tx2, tmp2::Tx22, v::Tx3, h::Tx4, hbar::Tx5, zerosvector::Tx6, u::Ty)
    LSMRDampenedAllocatedSolver{Tx0, Tx1, Tx2, Tx22, Tx3, Tx4, Tx5, Tx6, Ty}(preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, zerosvector, u)
end
function AbstractAllocatedSolver{Tx, Ty, Tf, TJ, Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::LevenbergMarquardt, solver::LSMR)
    preconditioner!, preconditioner = getpreconditioner(nls, optimizer, solver)
    LSMRDampenedAllocatedSolver(preconditioner!, preconditioner, _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x), _zeros(nls.x),  _zeros(nls.x), _zeros(nls.y))
end

function A_ldiv_B!(x, J, y, damp, A::LSMRDampenedAllocatedSolver)
    preconditioner!, preconditioner, tmp, tmp2, v, h, hbar, zerosvector, u =
            A.preconditioner!, A.preconditioner, A.tmp, A.tmp2, A.v, A.h, A.hbar, A.zerosvector, A.u
    # prepare x
    fill!(x, 0)

    # prepare b
    copy!(u, y)
    fill!(zerosvector, 0)
    b = DampenedVector(u, zerosvector)

    # prepare A
    fill!(tmp, 0)
    preconditioner!(x, J, damp, preconditioner)
    map!(sqrt, damp, damp)
    A = PreconditionedMatrix(DampenedMatrix(J, damp), preconditioner, tmp, tmp2)
    # solve
    x, ch = lsmr!(x, A, b, v, h, hbar, btol = 0.5)
    A_ldiv_B!(x, preconditioner, tmp)
    copy!(x, tmp)
    return x, ch.mvps
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests of formatters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Formatters" begin
    expected = """

    Col. 1
    Col. 2
    Col. 3
    Col. 4

    1
    false
    1
    1

    0
    true
    0
    0

    3
    false
    3
    3

    0
    true
    0
    0

    5
    false
    5
    5

    0
    true
    0
    0

"""

    formatter = (data,i,j) -> begin
        if j != 2
            return isodd(i) ? i : 0
        else
            return data
        end
    end

    result = pretty_table(String, data;
                          backend = Val(:html),
                          formatters = formatter,
                          standalone = false)

    @test result == expected
end