Rewrite IPM data structures (#62)

* rename directory

* Util function for generating vectors of given type

* Add an IPMData data structure and update IPM algorithm

* Add detailed timing within IPM optimization

* Rename files

* Conversion ProblemData --> IPMData and solution extraction

* Fix stopping criterion for Primal infeasibility

* Remove arithmetic for KKT solvers overview

* Update style convention in docstrings

* Update form of internal LP representation in docs

* Add toy example tutorial

* Bind IPMData to IPM optimizer and simplify dispatch rules

Author: mtanneau

Project.toml

@@ -14,12 +14,14 @@ QPSReader = "10f199a5-22af-520b-b891-7ce84a7b1bd0"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
 Krylov = "0.5.2"
 LDLFactorizations = "0.6"
 MathOptInterface = "0.9.5"
 QPSReader = "0.2"
+TimerOutputs = "0.5.6"
 julia = "1.3"
 
 [extras]

docs/Project.toml

@@ -1,3 +1,5 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Tulip = "6dd1b50a-3aae-11e9-10b5-ef983d2400fa"

docs/make.jl

@@ -6,6 +6,9 @@ makedocs(
     format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
     pages = [
         "Home" => "index.md",
+        "Tutorials" => Any[
+            "tutorials/lp_example.md",
+        ],
         "User manual" => Any[
             "Problem formulation" => "manual/formulation.md",
             "Solving linear systems" => "manual/linear_systems.md"

docs/src/manual/formulation.md

@@ -26,12 +26,11 @@ Internally, Tulip solves LPs of the form
     & c^{T} x + \ c_{0}\\
     s.t.
     & A x = b\\
-    & x \leq u\\
-    & x \geq 0
+    & l \leq x \leq u\\
     \end{array}
 ```
-where ``x, c, u \in \mathbb{R}^{n}``, ``A \in \mathbb{R}^{m \times n}`` and ``b \in \mathbb{R}^{m}``.
-Some ``u_{j}`` may may take infinite value, i.e., the corresponding variable ``x_{j}`` has no upper bound.
+where ``x, c \in \mathbb{R}^{n}``, ``A \in \mathbb{R}^{m \times n}``, ``b \in \mathbb{R}^{m}``,
+and  ``l, u \in (\mathbb{R} \cup \{-\infty, +\infty \})^{n}``, i.e., some bounds may be infinite.
 
 The original problem is automatically reformulated into standard form before the optimization is performed.
 This transformation is transparent to the user.

docs/src/manual/linear_systems.md

@@ -47,15 +47,15 @@ To enable the use of fast external libraries and/or specialized routines, the re
 
 Here is a list of currently supported linear solvers:
 
-| Linear solver type | `Tv` | System | Backend | Method |
-|:-------------------|:-----|:-------|:--------|:-------|
-| [`Dense_SymPosDef`](@ref) | `Real` | Normal equations | Dense / LAPACK | Cholesky
-| [`Cholmod_SymQuasDef`](@ref) | `Float64` | Augmented system | CHOLMOD | LDLᵀ
-| [`Cholmod_SymPosDef`](@ref) | `Float64` | Normal equations | CHOLMOD | Cholesky
-| [`LDLFact_SymQuasDef`](@ref) | `Real` | Augmented system | [LDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl) | LDLᵀ
-| [`KrylovSPDSolver`](@ref) | `Real` | Normal equations | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov
-| [`KrylovSIDSolver`](@ref) | `Real` | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov
-| [`KrylovSQDSolver`](@ref) | `Real` | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov
+| Linear solver type | System | Backend | Method |
+|:-------------------|:-------|:--------|:-------|
+| [`Dense_SymPosDef`](@ref) | Normal equations | Dense / LAPACK | Cholesky
+| [`Cholmod_SymQuasDef`](@ref) | Augmented system | CHOLMOD | LDLᵀ
+| [`Cholmod_SymPosDef`](@ref) | Normal equations | CHOLMOD | Cholesky
+| [`LDLFact_SymQuasDef`](@ref) | Augmented system | [LDLFactorizations.jl](https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl) | LDLᵀ
+| [`KrylovSPDSolver`](@ref) | Normal equations | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov
+| [`KrylovSIDSolver`](@ref) | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov
+| [`KrylovSQDSolver`](@ref) | Augmented system[^1] | [Krylov.jl](https://github.com/JuliaSmoothOptimizers/Krylov.jl) | Krylov
 
 [^1]: [`KrylovSIDSolver`](@ref)s view the augmented system as a symmetric indefinite system,
     while [`KrylovSQDSolver`](@ref)s exploit its 2x2 structure and quasi-definite property.

docs/src/reference/attributes.md

@@ -13,9 +13,9 @@ Attributes are queried using [`get_attribute`](@ref) and set using [`set_attribu
 | [`ModelName`](@ref)               | `String`  | Name of the model
 | [`NumberOfConstraints`](@ref)     | `Int`     | Number of constraints in the model
 | [`NumberOfVariables`](@ref)       | `Int`     | Number of variables in the model
-| [`ObjectiveValue`](@ref)          | `Tv`      | Objective value of the current primal solution
-| [`DualObjectiveValue`](@ref)      | `Tv`      | Objective value of the current dual solution
-| [`ObjectiveConstant`](@ref)       | `Tv`      | Value of the objective constant
+| [`ObjectiveValue`](@ref)          | `T`       | Objective value of the current primal solution
+| [`DualObjectiveValue`](@ref)      | `T`       | Objective value of the current dual solution
+| [`ObjectiveConstant`](@ref)       | `T`       | Value of the objective constant
 | [`ObjectiveSense`](@ref)          |           | Optimization sense
 | [`Status`](@ref)                  |           | Model status
 | [`BarrierIterations`](@ref)       | `Int`     | Number of barrier iterations
@@ -25,17 +25,17 @@ Attributes are queried using [`get_attribute`](@ref) and set using [`set_attribu
 
 | Name                               | Type     | Description                   
 |:-----------------------------------|:---------|:------------------------------
-| [`VariableLowerBound`](@ref)       | `Tv`     | Variable lower bound 
-| [`VariableUpperBound`](@ref)       | `Tv`     | Variable upper bound 
-| [`VariableObjectiveCoeff`](@ref)   | `Tv`     | Variable objective coefficient 
+| [`VariableLowerBound`](@ref)       | `T`      | Variable lower bound 
+| [`VariableUpperBound`](@ref)       | `T`      | Variable upper bound 
+| [`VariableObjectiveCoeff`](@ref)   | `T`      | Variable objective coefficient 
 | [`VariableName`](@ref)             | `String` | Variable name 
 
 ## Constraint attributes
 
 | Name                               | Type     | Description                   
 |:-----------------------------------|:---------|:------------------------------
-| [`ConstraintLowerBound`](@ref)     | `Tv`     | Constraint lower bound 
-| [`ConstraintUpperBound`](@ref)     | `Tv`     | Constraint upper bound
+| [`ConstraintLowerBound`](@ref)     | `T`      | Constraint lower bound 
+| [`ConstraintUpperBound`](@ref)     | `T`      | Constraint upper bound
 | [`ConstraintName`](@ref)           | `String` | Constraint name 
 
 
@@ -46,22 +46,22 @@ Attributes are queried using [`get_attribute`](@ref) and set using [`set_attribu
 
 ```@autodocs
 Modules = [Tulip]
-Pages = ["src/Core/attributes.jl"]
+Pages = ["src/attributes.jl"]
 Filter = t -> typeof(t) === DataType && t <: Tulip.AbstractModelAttribute
 ```
 
 ### Variable attributes
 
 ```@autodocs
 Modules = [Tulip]
-Pages = ["src/Core/attributes.jl"]
+Pages = ["src/attributes.jl"]
 Filter = t -> typeof(t) === DataType && t <: Tulip.AbstractVariableAttribute
 ```
 
 ### Constraint attributes
 
 ```@autodocs
 Modules = [Tulip]
-Pages = ["src/Core/attributes.jl"]
+Pages = ["src/attributes.jl"]
 Filter = t -> typeof(t) === DataType && t <: Tulip.AbstractConstraintAttribute
 ```

docs/src/tutorials/lp_example.md

@@ -0,0 +1,196 @@
+# Toy example
+
+Tulip can be accessed in 3 ways:
+through [JuMP](https://github.com/jump-dev/JuMP.jl),
+through [MathOptInterface](https://github.com/jump-dev/MathOptInterface.jl),
+or directly.
+
+This tutorial illustrates, for each case, how to build a model, solve it,
+and query the solution value.
+In all cases, we consider the small LP
+```math
+\begin{array}{rrrl}
+    (LP) \ \ \ 
+    \displaystyle Z^{*} = \min_{x, y} & -2x & - y\\
+    s.t.
+    &  x & - y & \geq -2\\
+    & 2x &-  y & \leq  4\\
+    &  x &+ 2y & \leq  7\\
+    &  x,&   y & \geq  0\\
+\end{array}
+```
+whose optimal value and solution are ``Z^{*} = -8`` and ``(x^{*}, y^{*}) = (3, 2)``.
+
+## JuMP
+
+```jldoctest; output = false
+using Printf
+using JuMP
+import Tulip
+
+# Instantiate JuMP model
+lp = Model(Tulip.Optimizer)
+
+# Create variables
+@variable(lp, x >= 0)
+@variable(lp, y >= 0)
+
+# Add constraints
+@constraint(lp, row1, x - y >= -2)
+@constraint(lp, row2, 2*x - y <= 4)
+@constraint(lp, row3, x + 2*y <= 7)
+
+# Set the objective
+@objective(lp, Min, -2*x - y)
+
+# Set some parameters
+set_optimizer_attribute(lp, "OutputLevel", 0)  # disable output
+set_optimizer_attribute(lp, "Presolve", 0)     # disable presolve
+
+# Solve the problem
+optimize!(lp)
+
+# Check termination status
+st = termination_status(lp)
+println("Termination status: $st")
+
+# Query solution value
+objval = objective_value(lp)
+x_ = value(x)
+y_ = value(y)
+
+@printf "Z* = %.4f\n" objval
+@printf "x* = %.4f\n" x_
+@printf "y* = %.4f\n" y_
+
+# output
+
+Termination status: OPTIMAL
+Z* = -8.0000
+x* = 3.0000
+y* = 2.0000
+```
+
+## MOI
+
+```jldoctest; output = false
+using Printf
+
+import MathOptInterface
+const MOI = MathOptInterface
+
+import Tulip
+
+lp = Tulip.Optimizer{Float64}()
+
+# Create variables
+x = MOI.add_variable(lp)
+y = MOI.add_variable(lp)
+
+# Set variable bounds
+MOI.add_constraint(lp, MOI.SingleVariable(x), MOI.GreaterThan(0.0))  # x >= 0
+MOI.add_constraint(lp, MOI.SingleVariable(y), MOI.GreaterThan(0.0))  # y >= 0
+
+# Add constraints
+row1 = MOI.add_constraint(lp,
+    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x, y]), 0.0),
+    MOI.GreaterThan(-2.0)
+)
+row2 = MOI.add_constraint(lp,
+    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, -1.0], [x, y]), 0.0),
+    MOI.LessThan(4.0)
+)
+row3 = MOI.add_constraint(lp,
+    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,  2.0], [x, y]), 0.0),
+    MOI.LessThan(7.0)
+) 
+
+# Set the objective
+MOI.set(lp,
+    MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
+    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-2.0, -1.0], [x, y]), 0.0)
+)
+MOI.set(lp, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+
+# Set some parameters
+MOI.set(lp, MOI.Silent(), true)               # disable output
+MOI.set(lp, MOI.RawParameter("Presolve"), 0)  # disable presolve
+
+# Solve the problem
+MOI.optimize!(lp)
+
+# Check status
+st = MOI.get(lp, MOI.TerminationStatus())
+println("Termination status: $st")
+
+# Query solution value
+objval = MOI.get(lp, MOI.ObjectiveValue())
+x_ = MOI.get(lp, MOI.VariablePrimal(), x)
+y_ = MOI.get(lp, MOI.VariablePrimal(), y)
+
+@printf "Z* = %.4f\n" objval
+@printf "x* = %.4f\n" x_
+@printf "y* = %.4f\n" y_
+
+# output
+
+Termination status: OPTIMAL
+Z* = -8.0000
+x* = 3.0000
+y* = 2.0000
+```
+
+## Tulip
+
+!!! warning
+    Tulip's low-level API should not be considered stable nor complete.
+    The recommended way to use Tulip is through JuMP/MOI as shown above.
+
+
+```jldoctest; output = false
+using Printf
+import Tulip
+
+# Instantiate Tulip object
+lp = Tulip.Model{Float64}()
+pb = lp.pbdata  # Internal problem data
+
+# Create variables
+x = Tulip.add_variable!(pb, Int[], Float64[], -2.0, 0.0, Inf, "x")
+y = Tulip.add_variable!(pb, Int[], Float64[], -1.0, 0.0, Inf, "y")
+
+# Add constraints
+row1 = Tulip.add_constraint!(pb, [x, y], [1.0, -1.0], -2.0, Inf, "row1")
+row2 = Tulip.add_constraint!(pb, [x, y], [2.0, -1.0], -Inf, 4.0, "row2")
+row3 = Tulip.add_constraint!(pb, [x, y], [1.0,  2.0], -Inf, 7.0, "row3")
+
+# Set the objective
+# Nothing to do here as objective is already declared
+
+# Set some parameters
+Tulip.set_parameter(lp, "OutputLevel", 0)  # disable output
+Tulip.set_parameter(lp, "Presolve", 0)     # disable presolve
+
+# Solve the problem
+Tulip.optimize!(lp)
+
+# Check termination status
+st = Tulip.get_attribute(lp, Tulip.Status())
+println("Termination status: $st")
+
+# Query solution value
+objval = Tulip.get_attribute(lp, Tulip.ObjectiveValue())
+x_ = lp.solution.x[x]
+y_ = lp.solution.x[y]
+
+@printf "Z* = %.4f\n" objval
+@printf "x* = %.4f\n" x_
+@printf "y* = %.4f\n" y_
+
+# output
+
+Termination status: Trm_Optimal
+Z* = -8.0000
+x* = 3.0000
+y* = 2.0000
+```