Add documentation for new features

docs/source/common_model_interface.md

@@ -39,11 +39,11 @@ continuous
 :return: the handle of the variable
 ```
 
-### Add multi-dimensional variables to the model as <project:#pyoptinterface.tupledict>
+### Add multidimensional variables to the model as <project:#pyoptinterface.tupledict>
 
 ```{py:function} model.add_variables(*coords, [lb=-inf, ub=+inf, domain=pyoptinterface.VariableDomain.Continuous, name=""])
 
-add a multi-dimensional variable to the model
+add a multidimensional variable to the model
 
 :param coords: the coordinates of the variable, can be a list of Iterables
 :param float lb: the lower bound of the variable, optional, defaults to $-\infty$
@@ -55,6 +55,22 @@ continuous
 :rtype: pyoptinterface.tupledict
 ```
 
+### Add multidimensional variables to the model as `numpy.ndarray`
+
+```{py:function} model.add_m_variables(shape, [lb=-inf, ub=+inf, domain=pyoptinterface.VariableDomain.Continuous, name=""])
+
+add a multidimensional variable to the model as `numpy.ndarray`
+
+:param shape: the shape of the variable, can be a tuple of integers or an integer
+:param float lb: the lower bound of the variable, optional, defaults to $-\infty$
+:param float ub: the upper bound of the variable, optional, defaults to $+\infty$
+:param pyoptinterface.VariableDomain domain: the domain of the variable, optional, defaults to 
+continuous
+:param str name: the name of the variable, optional
+:return: the multidimensional variable
+:rtype: numpy.ndarray
+```
+
 ### Get/set variable attributes
 
 ```{py:function} model.set_variable_attribute(var, attr, value)
@@ -137,6 +153,20 @@ pretty print an expression in a human-readable format
 - <project:#model.add_second_order_cone_constraint>
 - <project:#model.add_sos_constraint>
 
+### Add linear constraints as matrix form to the model
+
+```{py:function} model.add_m_linear_constraints(A, vars, sense, b, [name=""])
+
+add linear constraints as matrix form to the model $Ax \le b$ or $Ax = b$ or $Ax \ge b$
+
+:param A: the matrix of coefficients, can be a dense `numpy.ndarray` or a sparse matrix `scipy.sparse.sparray`
+:param vars: the variables in the constraints, can be a list or a 1-d `numpy.ndarray` returned by `add_m_variables`
+:param pyoptinterface.ConstraintSense sense: the sense of the constraints
+:param b: the right-hand side of the constraints, should be a 1-d `numpy.ndarray`
+:param str name: the name of the constraints, optional
+:return: the handles of linear constraints
+:rtype: numpy.ndarray
+```
 
 ### Get/set constraint attributes
 

docs/source/constraint.md

@@ -25,6 +25,16 @@ from pyoptinterface import copt
 model = copt.Model()
 ```
 
+## Constraint Sense
+
+The sense of a constraint can be one of the following:
+
+- `poi.Eq`: equal
+- `poi.Leq`: less than or equal
+- `poi.Geq`: greater than or equal
+
+They are the abbreviations of `poi.ConstraintSense.Equal`, `poi.ConstraintSense.LessEqual` or `poi.ConstraintSense.GreaterEqual` and can be used in the `sense` argument of the constraint creation functions.
+
 ## Linear Constraint
 It is defined as:
 
@@ -42,16 +52,15 @@ It can be added to the model using the `add_linear_constraint` method of the `Mo
 x = model.add_variable(name="x")
 y = model.add_variable(name="y")
 
-con = model.add_linear_constraint(2.0*x + 3.0*y, poi.ConstraintSense.LessEqual, 1.0)
+con = model.add_linear_constraint(2.0*x + 3.0*y, poi.Leq, 1.0)
 ```
 
 ```{py:function} model.add_linear_constraint(expr, sense, rhs, [name=""])
 
 add a linear constraint to the model
 
 :param expr: the expression of the constraint
-:param pyoptinterface.ConstraintSense sense: the sense 
-of the constraint, which can be `GreaterEqual`, `Equal`, or `LessEqual`
+:param pyoptinterface.ConstraintSense sense: the sense of the constraint
 :param float rhs: the right-hand side of the constraint
 :param str name: the name of the constraint, optional
 :return: the handle of the constraint
@@ -88,8 +97,7 @@ con = model.add_quadratic_constraint(expr, poi.ConstraintSense.LessEqual, 1.0)
 add a quadratic constraint to the model
 
 :param expr: the expression of the constraint
-:param pyoptinterface.ConstraintSense sense: the sense 
-of the constraint, which can be `GreaterEqual`, `Equal`, or `LessEqual`
+:param pyoptinterface.ConstraintSense sense: the sense of the constraint, which can be `GreaterEqual`, `Equal`, or `LessEqual`
 :param float rhs: the right-hand side of the constraint
 :param str name: the name of the constraint, optional
 :return: the handle of the constraint
@@ -193,6 +201,8 @@ standard [constraint attributes](#pyoptinterface.ConstraintAttribute):
     - float
 *   - Dual
     - float
+*   - IIS
+    - bool
 :::
 
 The most common attribute we will use is the `Dual` attribute, which represents the dual multiplier of the constraint after optimization.
@@ -235,3 +245,26 @@ model.set_normalized_rhs(con, 2.0)
 # modify the coefficient of the linear part of the constraint
 model.set_normalized_coefficient(con, x, 2.0)
 ```
+
+## Create constraint with comparison operator
+
+In other modeling languages, we can create a constraint with a comparison operator, like:
+
+```python
+model.addConsr(x + y <= 1)
+```
+
+This is quite convenient, so PyOptInterface now supports to create constraint with comparison operators `<=`, `==`, `>=` as a shortcut to create a linear or quadratic constraint.
+
+```{code-cell}
+model.add_linear_constraint(x + y <= 1)
+model.add_linear_constraint(x <= y)
+model.add_quadratic_constraint(x*x + y*y <= 1)
+```
+
+:::{note}
+
+Creating constraint with comparison operator may cause performance issue especially the left-hand side and right-hand side of the constraint are complex expressions. PyOptInterface needs to create a new expression by subtracting the right-hand side from the left-hand side, which may be time-consuming.
+
+If that becomes the bottleneck of performance, it is recommended to construct the left-hand side expression with `ExprBuilder` and call `add_linear_constraint` or `add_quadratic_constraint` method to create constraints explicitly.
+:::

docs/source/faq.md

@@ -34,110 +34,4 @@ model = gurobi.Model(env)
 
 In YALMIP, you can use the matrix form $Ax \leq b$ to add linear constraints, which is quite convenient.
 
-In PyOptInterface, you can use the following code to add linear constraints in matrix form:
-
-```python
-import pyoptinterface as poi
-from pyoptinterface import gurobi
-
-import numpy as np
-from scipy.sparse import csr_array, sparray, eye_array
-
-
-def iterate_sparse_matrix_rows(A):
-    """
-    Iterate over rows of a sparse matrix and get non-zero elements for each row.
-
-    A is a 2-dimensional scipy sparse matrix
-    isinstance(A, scipy.sparse.sparray) = True and A.ndim = 2
-    """
-    if not isinstance(A, csr_array):
-        A = csr_array(A)  # Convert to CSR format if not already
-
-    for i in range(A.shape[0]):
-        row_start = A.indptr[i]
-        row_end = A.indptr[i + 1]
-        row_indices = A.indices[row_start:row_end]
-        row_data = A.data[row_start:row_end]
-        yield row_indices, row_data
-
-
-def add_matrix_constraints(model, A, x, sense, b):
-    """
-    add constraints Ax <= / = / >= b
-
-    A is a 2-dimensional numpy array or scipy sparse matrix
-    x is an iterable of variables
-    sense is one of (poi.Leq, poi.Eq, poi.Geq)
-    b is an iterable of values or a single scalar
-    """
-
-    is_ndarray = isinstance(A, np.ndarray)
-    is_sparse = isinstance(A, sparray)
-
-    if not is_ndarray and not is_sparse:
-        raise ValueError("A must be a numpy array or scipy.sparse array")
-
-    ndim = A.ndim
-    if ndim != 2:
-        raise ValueError("A must be a 2-dimensional array")
-
-    M, N = A.shape
-
-    # turn x into a list if x is an iterable
-    if isinstance(x, poi.tupledict):
-        x = x.values()
-    x = list(x)
-
-    if len(x) != N:
-        raise ValueError("x must have length equal to the number of columns of A")
-
-    # check b
-    if np.isscalar(b):
-        b = np.full(M, b)
-    elif len(b) != M:
-        raise ValueError("b must have length equal to the number of rows of A")
-
-    constraints = []
-
-    if is_ndarray:
-        for i in range(M):
-            expr = poi.ScalarAffineFunction()
-            row = A[i]
-            for coef, var in zip(row, x):
-                expr.add_term(var, coef)
-            con = model.add_linear_constraint(expr, sense, b[i])
-            constraints.append(con)
-    elif is_sparse:
-        for (row_indices, row_data), rhs in zip(iterate_sparse_matrix_rows(A), b):
-            expr = poi.ScalarAffineFunction()
-            for j, coef in zip(row_indices, row_data):
-                expr.add_term(x[j], coef)
-            con = model.add_linear_constraint(expr, sense, rhs)
-            constraints.append(con)
-
-    return constraints
-
-
-def main():
-    model = gurobi.Model()
-    N = 200
-    x = model.add_variables(range(N))
-    A = np.eye(N)
-    ub = 3.0
-    lb = 1.0
-    A_sparse = eye_array(N)
-    add_matrix_constraints(model, A, x, poi.Leq, ub)
-    add_matrix_constraints(model, A_sparse, x, poi.Geq, lb)
-
-    obj = poi.quicksum(x)
-    model.set_objective(obj)
-    model.optimize()
-
-    obj_value = model.get_model_attribute(poi.ModelAttribute.ObjectiveValue)
-    print("Objective value: ", obj_value)
-
-
-if __name__ == "__main__":
-    main()
-```
+In PyOptInterface, you can use [`model.add_m_linear_constraints`](<project:#model.add_m_linear_constraints>) to add linear constraints in matrix form.

docs/source/getting_started.md

@@ -245,7 +245,7 @@ con = model.add_linear_constraint(x1+x2, poi.ConstraintSense.Equal, 1, name="con
 ```
 `model.add_linear_constraint` adds a linear constraint to the model.
 - The first argument `x1+x2` is the left-hand side of the constraint.
-- The second argument is the sense of the constraint. It can be `poi.ConstraintSense.Equal`, `poi.ConstraintSense.LessEqual` or `poi.ConstraintSense.GreaterEqual`.
+- The second argument is the sense of the constraint. It can be `poi.ConstraintSense.Equal`, `poi.ConstraintSense.LessEqual` or `poi.ConstraintSense.GreaterEqual` which can also be written as `poi.Eq`, `poi.Leq`, and `poi.Geq`.
 - The third argument is the right-hand side of the constraint. It must be a constant.
 - The fourth argument is optional and can be used to specify the name of the constraint.
 

docs/source/index.md

@@ -20,6 +20,7 @@ container.md
 numpy.md
 structure.md
 common_model_interface.md
+infeasibility.md
 callback.md
 gurobi.md
 copt.md

docs/source/infeasibility.md

@@ -0,0 +1,68 @@
+---
+file_format: mystnb
+kernelspec:
+  name: python3
+---
+
+# Infeasibility Analysis
+
+The optimization model is not ways feasible, and the optimizer may tell us some information about the infeasibility to diagnose the problem. There are two ways to  handle the infeasibilities:
+
+- Find the IIS (Irreducible Infeasible Set) to identify the minimal set of constraints that cause the infeasibility.
+- Relax the constraints and solve a weaker problem to find out which constraints are violated and how much.
+
+PyOptInterface currently supports the first method to find the IIS (only with Gurobi and COPT). The following code snippet shows how to find the IIS of an infeasible model:
+
+```{code-cell}
+import pyoptinterface as poi
+from pyoptinterface import copt
+
+model = copt.Model()
+
+x = model.add_variable(lb=0.0, name="x")
+y = model.add_variable(lb=0.0, name="y")
+
+con1 = model.add_linear_constraint(x + y, poi.Geq, 5.0)
+con2 = model.add_linear_constraint(x + 2 * y, poi.Leq, 1.0)
+
+model.set_objective(x)
+
+model.computeIIS()
+
+con1_iis = model.get_constraint_attribute(con1, poi.ConstraintAttribute.IIS)
+con2_iis = model.get_constraint_attribute(con2, poi.ConstraintAttribute.IIS)
+
+print(f"Constraint 1 IIS: {con1_iis}")
+print(f"Constraint 2 IIS: {con2_iis}")
+```
+
+This code snippet creates an infeasible model with two constraints and finds the IIS of the model. Obviously, the constraints are contradictory because `x + 2 * y <= 1` and `x + y >= 5` cannot be satisfied at the same time when `x` and `y` are non-negative. The optimizer will detect that the model is infeasible and return the IIS, which is the set of constraints that cause the infeasibility. We can query whether a constraint is in the IIS by calling `get_constraint_attribute` with the `ConstraintAttribute.IIS` attribute.
+
+Sometimes, the bounds of the variables are not consistent with the constraints, and we need to query the IIS of the bounds of variables by calling `get_variable_attribute` with the `VariableAttribute.IISLowerBound` and `VariableAttribute.IISUpperBound` attributes.
+
+The following code snippet shows how to tell if the bounds of a variable are in the IIS:
+
+```{code-cell}
+model = copt.Model()
+
+x = model.add_variable(lb=0.0, ub=2.0, name="x")
+y = model.add_variable(lb=0.0, ub=3.0, name="y")
+
+con1 = model.add_linear_constraint(x + y, poi.Geq, 6.0)
+
+model.set_objective(x)
+
+model.computeIIS()
+
+con1_iis = model.get_constraint_attribute(con1, poi.ConstraintAttribute.IIS)
+x_lb_iis = model.get_variable_attribute(x, poi.VariableAttribute.IISLowerBound)
+x_ub_iis = model.get_variable_attribute(x, poi.VariableAttribute.IISUpperBound)
+y_lb_iis = model.get_variable_attribute(y, poi.VariableAttribute.IISLowerBound)
+y_ub_iis = model.get_variable_attribute(y, poi.VariableAttribute.IISUpperBound)
+
+print(f"Constraint 1 IIS: {con1_iis}")
+print(f"Variable x lower bound IIS: {x_lb_iis}")
+print(f"Variable x upper bound IIS: {x_ub_iis}")
+print(f"Variable y lower bound IIS: {y_lb_iis}")
+print(f"Variable y upper bound IIS: {y_ub_iis}")
+```

docs/source/numpy.md

@@ -4,13 +4,15 @@ kernelspec:
   name: python3
 ---
 
-# Numpy Container and N-queens Problem
+# Matrix Modeling
 
-In the previous [container](container.md) section, we have introduced the `tupledict` container to store and manipulate multi-dimensional data.
+In the previous [container](container.md) section, we have introduced the `tupledict` container to store and manipulate multidimensional data.
 
-However, due to the Bring Your Own Container (BYOC) principle, variables and constraints in PyOptInterface can just simple Python objects that can be stored in Numpy `ndarrays` directly as a multi-dimensional array, and you can enjoy the features of Numpy such like [fancy-indexing](https://numpy.org/doc/stable/user/basics.indexing.html) automatically. 
+However, due to the Bring Your Own Container (BYOC) principle, variables and constraints in PyOptInterface can just simple Python objects that can be stored in Numpy `ndarray` directly as a multidimensional array, and you can enjoy the features of Numpy such like [fancy-indexing](https://numpy.org/doc/stable/user/basics.indexing.html) automatically.
 
-We will use N-queens problem as example to show how to use Numpy `ndarrays` as container to store 2-dimensional variables and construct optimization model.
+## N-queen problem
+
+We will use N-queens problem as example to show how to use Numpy `ndarray` as container to store 2-dimensional variables and construct optimization model.
 
 Firstly, we import the necessary modules:
 
@@ -63,3 +65,32 @@ x_value = get_v(x)
 
 print(x_value.astype(int))
 ```
+
+## Built-in functions to add variables and constraints as Numpy `ndarray`
+
+Although you can construct the `ndarray` of variables and constraints manually, PyOptInterface provides built-in functions to simplify the process. The following code snippet shows how to use the built-in functions to add variables and constraints as Numpy `ndarray`:
+
+```{code-cell}
+model = highs.Model()
+
+x = model.add_m_variables(N)
+
+A = np.eye(N)
+b_ub = np.ones(N)
+b_lb = np.ones(N)
+
+model.add_m_linear_constraints(A, x, poi.Leq, b_ub)
+model.add_m_linear_constraints(A, x, poi.Geq, b_lb)
+
+model.set_objective(poi.quicksum(x))
+
+model.optimize()
+```
+
+Here we use two built-in functions `add_m_variables` and `add_m_linear_constraints` to add variables and constraints as Numpy `ndarray` respectively.
+
+The reference of these functions are listed in <project:#model.add_m_variables> and <project:#model.add_m_linear_constraints>.
+
+`add_m_variables` returns a `ndarray` of variables with the specified shape.
+
+`add_m_linear_constraints` adds multiple linear constraints to the model at once formulated as $Ax \le b$ or $Ax = b$ or $Ax \ge b$ where the matrix $A$ can be a dense `numpy.ndarray` or a sparse matrix `scipy.sparse.sparray`.