Repository: Qiskit/qiskit-optimization (python)

Issue:
Special patterns of inequality to penalty converters
<!-- ⚠️ If you do not respect this template, your issue will be closed -->
<!-- ⚠️ Make sure to browse the opened and closed issues to confirm this idea does not exist. -->

### What is the expected enhancement?

When users solve an optimization problem with Qiskit Aqua, they apply `InequalityToEquality` to convert inequality constraints into equality constraints by introducing slack variables, and then apply `LinearEqualityToPenalty` to translate the constraints into penalties of the objective function of QUBO. But, there are some special patterns are known that does not require slack variables. The objective of this issue is to implement such special converters.

[A Tutorial on Formulating and Using QUBO Models](https://arxiv.org/abs/1811.11538) introduces examples of the special patterns in page 10 as follows.
![image](https://user-images.githubusercontent.com/31178928/92363149-08fee900-f12c-11ea-9ca7-74b1a4b5f26d.png)



Pull request:
Special patterns of inequality constraint to penalty converters
<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add it in the CHANGELOG file under Unreleased section.
⚠️ If your pull request fixes an open issue, please link to the issue.

✅ I have added the tests to cover my changes.
✅ I have updated the documentation accordingly.
✅ I have read the CONTRIBUTING document.
-->

### Summary
Resolves #40. For details about previous discussion, see #98 

### Details and comments
Owned by Amano-san and me

[A Tutorial on Formulating and Using QUBO Models](https://arxiv.org/abs/1811.11538) introduces examples of the special patterns in page 10 as follows.
![image](https://user-images.githubusercontent.com/31178928/92363149-08fee900-f12c-11ea-9ca7-74b1a4b5f26d.png)


Patch:
diff --git a/.pylintdict b/.pylintdict
index b613937..cdd2c65 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -57,11 +57,13 @@ farhi
 fmin
 formatter
 func
+fred
 fval
 fx
 gambella
 geq
 getter
+glover
 goemans
 goldstone
 grover
@@ -155,6 +157,7 @@ representable
 rhobeg
 rhoend
 rhs
+rightarrow
 ronagh
 rtype
 runtime
diff --git a/qiskit_optimization/converters/__init__.py b/qiskit_optimization/converters/__init__.py
index f796e27..7b56e92 100644
--- a/qiskit_optimization/converters/__init__.py
+++ b/qiskit_optimization/converters/__init__.py
@@ -40,6 +40,7 @@ Converters
    InequalityToEquality
    IntegerToBinary
    LinearEqualityToPenalty
+   LinearInequalityToPenalty
    MaximizeToMinimize
    MinimizeToMaximize
    QuadraticProgramToQubo
@@ -49,6 +50,7 @@ Converters
 from .integer_to_binary import IntegerToBinary
 from .inequality_to_equality import InequalityToEquality
 from .linear_equality_to_penalty import LinearEqualityToPenalty
+from .linear_inequality_to_penalty import LinearInequalityToPenalty
 from .flip_problem_sense import MaximizeToMinimize
 from .flip_problem_sense import MinimizeToMaximize
 from .quadratic_program_to_qubo import QuadraticProgramToQubo
@@ -58,6 +60,7 @@ __all__ = [
     "InequalityToEquality",
     "IntegerToBinary",
     "LinearEqualityToPenalty",
+    "LinearInequalityToPenalty",
     "MaximizeToMinimize",
     "MinimizeToMaximize",
     "QuadraticProgramConverter",
diff --git a/qiskit_optimization/converters/linear_inequality_to_penalty.py b/qiskit_optimization/converters/linear_inequality_to_penalty.py
new file mode 100644
index 0000000..0bf08ab
--- /dev/null
+++ b/qiskit_optimization/converters/linear_inequality_to_penalty.py
@@ -0,0 +1,356 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2020, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Converter to convert a problem with inequality constraints to unconstrained with penalty terms."""
+
+import logging
+from typing import Optional, Union, Tuple, List, Dict
+
+import numpy as np
+
+from .quadratic_program_converter import QuadraticProgramConverter
+from ..exceptions import QiskitOptimizationError
+from ..problems.constraint import Constraint, ConstraintSense
+from ..problems.quadratic_objective import QuadraticObjective
+from ..problems.quadratic_program import QuadraticProgram
+from ..problems.variable import Variable
+
+logger = logging.getLogger(__name__)
+
+
+class LinearInequalityToPenalty(QuadraticProgramConverter):
+    r"""Convert linear inequality constraints to penalty terms of the objective function.
+
+    There are some linear constraints which do not require slack variables to
+    construct penalty terms [1]. This class supports the following inequality constraints.
+
+    .. math::
+
+        \begin{array}{}
+        \text { Inequality constraint } & & \text { Penalty term } \\
+        x \leq y & \rightarrow  & P(x-x y) \\
+        x \geq y & \rightarrow  & P(y-x y) \\
+        \sum_{i=1}^n x_i \leq 1, n \geq 2 & \rightarrow & P \sum_{i, j : i < j} x_i x_j\\
+        \sum_{i=1}^n x_i \geq n-1, n \geq 2 & \rightarrow & P \sum_{i, j : i < j} (1 - x_i) (1 - x_j)
+        \end{array}
+
+    Note that x, y, z and :math:`x_i` are binary variables, and P is a penalty factor,
+    where the value of P is automatically determined or supplied by users.
+
+    If constraints match with any of the patterns, they are converted into penalty terms and added
+    to the objective function. Otherwise, constraints are kept as is.
+
+    References:
+        [1]: Fred Glover, et al. (2019),
+             A Tutorial on Formulating and Using QUBO Models,
+             `arXiv:1811.11538 <https://arxiv.org/abs/1811.11538>`_.
+
+    """
+
+    def __init__(self, penalty: Optional[float] = None) -> None:
+        """
+        Args:
+            penalty: Penalty factor to scale equality constraints that are added to objective.
+                     If None is passed, a penalty factor will be automatically calculated on
+                     every conversion.
+        """
+
+        self._src_num_vars: Optional[int] = None
+        self._dst: Optional[QuadraticProgram] = None
+        self._penalty: Optional[float] = penalty
+        self._should_define_penalty: bool = penalty is None
+
+    def convert(self, problem: QuadraticProgram) -> QuadraticProgram:
+        r"""Convert inequality constraints into penalty terms of the objective function.
+
+        This methods converts the following patterns where x, y, and :math:`x_i` are binary variables
+        and P is a penalty factor.
+
+        .. math::
+
+            \begin{array}{}
+            \text { Inequality constraint } & & \text { Penalty term } \\
+            x \leq y & \rightarrow  & P(x-x y) \\
+            x \geq y & \rightarrow  & P(y-x y) \\
+            \sum_{i=1}^n x_i \leq 1, n \geq 2 & \rightarrow & P \sum_{i, j : i < j} x_i x_j\\
+            \sum_{i=1}^n x_i \geq n-1, n \geq 2 & \rightarrow & P \sum_{i, j : i < j} (1 - x_i) (1 - x_j)
+            \end{array}
+
+        Args:
+            problem: The problem to be solved.
+
+        Returns:
+            The converted problem
+
+        Raises:
+            QiskitOptimizationError: If an unsupported-type variable exists.
+        """
+
+        # create empty QuadraticProgram model
+        self._src_num_vars = problem.get_num_vars()
+        self._dst = QuadraticProgram(name=problem.name)
+
+        # If no penalty was given, set the penalty coefficient by _auto_define_penalty()
+        if self._should_define_penalty:
+            penalty = self._auto_define_penalty(problem)
+        else:
+            penalty = self._penalty
+
+        # Set variables
+        for x in problem.variables:
+            if x.vartype == Variable.Type.CONTINUOUS:
+                self._dst.continuous_var(x.lowerbound, x.upperbound, x.name)
+            elif x.vartype == Variable.Type.BINARY:
+                self._dst.binary_var(x.name)
+            elif x.vartype == Variable.Type.INTEGER:
+                self._dst.integer_var(x.lowerbound, x.upperbound, x.name)
+            else:
+                raise QiskitOptimizationError(f"Unsupported vartype: {x.vartype}")
+
+        # get original objective terms
+        offset = problem.objective.constant
+        linear = problem.objective.linear.to_dict()
+        quadratic = problem.objective.quadratic.to_dict()
+        sense = problem.objective.sense.value
+
+        # convert linear constraints into penalty terms
+        for constraint in problem.linear_constraints:
+
+            # special constraint check function here
+            if not self._is_matched_constraint(problem, constraint):
+                self._dst.linear_constraint(
+                    constraint.linear.coefficients,
+                    constraint.sense,
+                    constraint.rhs,
+                    constraint.name,
+                )
+                continue
+
+            conv_offset, conv_linear, conv_quadratic, varmap = self._conversion_table(constraint)
+
+            # constant part
+            offset += sense * penalty * conv_offset
+
+            # linear parts of penalty
+            for j, j_2 in varmap.items():
+                # if j already exists in the linear terms dic, add a penalty term
+                # into existing value else create new key and value in the linear_term dict
+
+                if conv_linear[j] != 0:
+                    linear[j_2] = linear.get(j_2, 0.0) + sense * penalty * conv_linear[j]
+
+            # quadratic parts of penalty
+            for j, j_2 in varmap.items():
+                for k in range(j, len(varmap)):
+                    # if j and k already exist in the quadratic terms dict,
+                    # add a penalty term into existing value
+                    # else create new key and value in the quadratic term dict
+
+                    if conv_quadratic[j][k] != 0:
+                        tup = (j_2, varmap[k])
+                        quadratic[tup] = (
+                            quadratic.get(tup, 0.0) + sense * penalty * conv_quadratic[j][k]
+                        )
+
+        # Copy quadratic_constraints
+        for quadratic_constraint in problem.quadratic_constraints:
+            self._dst.quadratic_constraint(
+                quadratic_constraint.linear.coefficients,
+                quadratic_constraint.quadratic.coefficients,
+                quadratic_constraint.sense,
+                quadratic_constraint.rhs,
+                quadratic_constraint.name,
+            )
+
+        if problem.objective.sense == QuadraticObjective.Sense.MINIMIZE:
+            self._dst.minimize(offset, linear, quadratic)
+        else:
+            self._dst.maximize(offset, linear, quadratic)
+
+        # Update the penalty to the one just used
+        self._penalty = penalty
+
+        return self._dst
+
+    @staticmethod
+    def _conversion_table(
+        constraint,
+    ) -> Tuple[int, np.ndarray, np.ndarray, Dict[int, int]]:
+        """Construct conversion matrix for special constraint.
+
+        Returns:
+            Return conversion table which is used to construct
+            penalty term in main function.
+
+        Raises:
+            QiskitOptimizationError: if the constraint is invalid.
+        """
+
+        vars_dict = constraint.linear.to_dict()
+        coeffs = list(vars_dict.values())
+        varmap = dict(enumerate(vars_dict.keys()))
+        rhs = constraint.rhs
+        sense = constraint.sense
+
+        num_vars = len(vars_dict)
+
+        # initialize return values, these are used for converted offset, linear
+        # and quadratic terms
+        offset = 0
+        linear = np.zeros(num_vars, dtype=int)
+        quadratic = np.zeros((num_vars, num_vars), dtype=int)
+
+        # rhs = num_vars - 1 correspond to multiple variable with >= n - 1 case.
+        if sense == ConstraintSense.GE and rhs == num_vars - 1:
+            # x_1 + ... + x_n >= n - 1
+            # The number of offset is combination ( nC2 )
+            offset = num_vars * (num_vars - 1) // 2
+            linear = np.full(num_vars, 1 - num_vars, dtype=int)
+            quadratic = np.triu(np.ones((num_vars, num_vars), dtype=int), k=1)
+        elif sense == ConstraintSense.LE and rhs == 1:
+            # x_1 + ... + x_n <= 1
+            quadratic = np.triu(np.ones((num_vars, num_vars), dtype=int), k=1)
+        elif rhs == 0:
+            if num_vars != 2:
+                raise QiskitOptimizationError(
+                    f"Internal error: invalid number of variables {num_vars} {constraint.name}"
+                )
+            quadratic = np.array([[0, -1], [0, 0]])
+            if sense == ConstraintSense.GE:
+                # x >= y case
+                if coeffs[0] < 0.0:
+                    linear[0] = 1
+                else:
+                    linear[1] = 1
+            elif sense == ConstraintSense.LE:
+                # x <= y case
+                if coeffs[0] > 0.0:
+                    linear[0] = 1
+                else:
+                    linear[1] = 1
+        else:
+            raise QiskitOptimizationError(f"Internal error: invalid constraint {constraint.name}")
+
+        return offset, linear, quadratic, varmap
+
+    @staticmethod
+    def _is_matched_constraint(problem, constraint) -> bool:
+        """Determine if constraint is special or not.
+
+        Returns:
+            True: when constraint is special
+            False: when constraint is not special
+        """
+
+        params = constraint.linear.to_dict()
+        num_vars = len(params)
+        rhs = constraint.rhs
+        sense = constraint.sense
+        coeff_array = np.array(list(params.values()))
+
+        # Binary parameter?
+        if any(problem.variables[i].vartype != Variable.Type.BINARY for i in params.keys()):
+            return False
+
+        if num_vars == 2 and rhs == 0:
+            if sense in (Constraint.Sense.LE, Constraint.Sense.GE):
+                # x-y<=0
+                # x-y>=0
+                return coeff_array.min() == -1.0 and coeff_array.max() == 1.0
+        elif num_vars >= 2:
+            if sense == Constraint.Sense.LE and rhs == 1:
+                if all(i == 1 for i in params.values()):
+                    # x1+x2+...<=1
+                    return True
+            elif sense == Constraint.Sense.GE and rhs == num_vars - 1:
+                if all(i == 1 for i in params.values()):
+                    # x1+x2+...>=n-1
+                    return True
+
+        return False
+
+    @staticmethod
+    def _auto_define_penalty(problem) -> float:
+        """Automatically define the penalty coefficient.
+
+        Returns:
+            Return the minimum valid penalty factor calculated
+            from the upper bound and the lower bound of the objective function.
+            If a constraint has a float coefficient,
+            return the default value for the penalty factor.
+        """
+
+        default_penalty = 1e5
+
+        # Check coefficients of constraints.
+        # If a constraint has a float coefficient, return the default value for the penalty factor.
+        terms = []
+        for constraint in problem.linear_constraints:
+            terms.append(constraint.rhs)
+            terms.extend(constraint.linear.to_array().tolist())
+        if any(isinstance(term, float) and not term.is_integer() for term in terms):
+            logger.warning(
+                "Warning: Using %f for the penalty coefficient because "
+                "a float coefficient exists in constraints. \n"
+                "The value could be too small. "
+                "If so, set the penalty coefficient manually.",
+                default_penalty,
+            )
+            return default_penalty
+
+        lin_b = problem.objective.linear.bounds
+        quad_b = problem.objective.quadratic.bounds
+        return 1.0 + (lin_b.upperbound - lin_b.lowerbound) + (quad_b.upperbound - quad_b.lowerbound)
+
+    def interpret(self, x: Union[np.ndarray, List[float]]) -> np.ndarray:
+        """Convert the result of the converted problem back to that of the original problem
+
+        Args:
+            x: The result of the converted problem or the given result in case of FAILURE.
+
+        Returns:
+            The result of the original problem.
+
+        Raises:
+            QiskitOptimizationError: if the number of variables in the result differs from
+                                     that of the original problem.
+        """
+
+        if len(x) != self._src_num_vars:
+            raise QiskitOptimizationError(
+                f"The number of variables in the passed result ({len(x)}) differs from "
+                f"that of the original problem ({self._src_num_vars})."
+            )
+        return np.asarray(x)
+
+    @property
+    def penalty(self) -> Optional[float]:
+        """Returns the penalty factor used in conversion.
+
+        Returns:
+            The penalty factor used in conversion.
+        """
+
+        return self._penalty
+
+    @penalty.setter
+    def penalty(self, penalty: Optional[float]) -> None:
+        """Set a new penalty factor.
+
+        Args:
+            penalty: The new penalty factor.
+                     If None is passed, a penalty factor will be automatically calculated
+                     on every conversion.
+        """
+
+        self._penalty = penalty
+        self._should_define_penalty = penalty is None
diff --git a/qiskit_optimization/converters/quadratic_program_to_qubo.py b/qiskit_optimization/converters/quadratic_program_to_qubo.py
index 424bbda..3ca6d48 100644
--- a/qiskit_optimization/converters/quadratic_program_to_qubo.py
+++ b/qiskit_optimization/converters/quadratic_program_to_qubo.py
@@ -20,6 +20,7 @@ from ..converters.flip_problem_sense import MaximizeToMinimize
 from ..converters.inequality_to_equality import InequalityToEquality
 from ..converters.integer_to_binary import IntegerToBinary
 from ..converters.linear_equality_to_penalty import LinearEqualityToPenalty
+from ..converters.linear_inequality_to_penalty import LinearInequalityToPenalty
 from ..exceptions import QiskitOptimizationError
 from ..problems.quadratic_program import QuadraticProgram
 from .quadratic_program_converter import QuadraticProgramConverter
@@ -44,9 +45,10 @@ class QuadraticProgramToQubo(QuadraticProgramConverter):
                 If None is passed, a penalty factor will be automatically calculated on every
                 conversion.
         """
-
         self._penalize_lin_eq_constraints = LinearEqualityToPenalty(penalty=penalty)
+        self._penalize_lin_ineq_constraints = LinearInequalityToPenalty(penalty=penalty)
         self._converters = [
+            self._penalize_lin_ineq_constraints,
             InequalityToEquality(mode="integer"),
             IntegerToBinary(),
             self._penalize_lin_eq_constraints,
@@ -160,4 +162,5 @@ class QuadraticProgramToQubo(QuadraticProgramConverter):
                      If None is passed, a penalty factor will be automatically calculated on every
                      conversion.
         """
+        self._penalize_lin_ineq_constraints.penalty = penalty
         self._penalize_lin_eq_constraints.penalty = penalty
diff --git a/releasenotes/notes/linear-inequality-to-penalty-70b5277b8edda012.yaml b/releasenotes/notes/linear-inequality-to-penalty-70b5277b8edda012.yaml
new file mode 100644
index 0000000..5d012d8
--- /dev/null
+++ b/releasenotes/notes/linear-inequality-to-penalty-70b5277b8edda012.yaml
@@ -0,0 +1,16 @@
+---
+features:
+  - |
+    Introduced a new converter class :class:`qiskit_optimization.converters.LinearInequalityToPenalty`.
+    It converts the following inequality constraints to penalty terms where x, y, :math:`x_i`
+    are binary variables and P is a penalty factor.
+
+    .. math::
+
+        \begin{array}{}
+        \text { Inequality constraint } & & \text { Penalty term } \\
+        x \leq y & \rightarrow  & P(x-x y) \\
+        x \geq y & \rightarrow  & P(y-x y) \\
+        \sum_{i=1}^n x_i \leq 1, n \geq 2 & \rightarrow & P \sum_{i, j : i < j} x_i x_j\\
+        \sum_{i=1}^n x_i \geq n-1, n \geq 2 & \rightarrow & P \sum_{i, j : i < j} (1 - x_i) (1 - x_j)
+        \end{array}


Test patch:
diff --git a/test/converters/test_linear_inequality_to_penalty.py b/test/converters/test_linear_inequality_to_penalty.py
new file mode 100644
index 0000000..a394d0f
--- /dev/null
+++ b/test/converters/test_linear_inequality_to_penalty.py
@@ -0,0 +1,566 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+""" Test LinearInequalityToPenalty converter """
+
+import unittest
+from test.optimization_test_case import QiskitOptimizationTestCase
+
+from qiskit_optimization import QuadraticProgram
+from qiskit_optimization.converters import LinearInequalityToPenalty, QuadraticProgramToQubo
+from qiskit_optimization.problems import Constraint
+
+
+class TestLinearInequalityToPenalty(QiskitOptimizationTestCase):
+    """Test Converters"""
+
+    def test_linear_inequality_to_penalty1(self):
+        """Test special constraint to penalty x+y <= 1 -> P(x*y)"""
+
+        op = QuadraticProgram()
+        lip = LinearInequalityToPenalty()
+
+        op.binary_var(name="x")
+        op.binary_var(name="y")
+        # Linear constraints
+        linear_constraint = {"x": 1, "y": 1}
+        op.linear_constraint(linear_constraint, Constraint.Sense.LE, 1, "P(xy)")
+
+        # Test with no max/min
+        with self.subTest("No max/min"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            lip.penalty = 1
+            quadratic = {("x", "y"): lip.penalty}
+            op2 = lip.convert(op)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test maximize
+        with self.subTest("Maximize"):
+            linear = {"x": 2, "y": 1}
+            op.maximize(linear=linear)
+            lip.penalty = 5
+            quadratic = {("x", "y"): -1 * lip.penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test minimize
+        with self.subTest("Minimize"):
+            linear = {"x": 2, "y": 1}
+            op.minimize(linear=linear)
+            lip.penalty = 5
+            quadratic = {("x", "y"): lip.penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test combination
+        with self.subTest("Combination"):
+            op = QuadraticProgram()
+            lip = LinearInequalityToPenalty()
+
+            op.binary_var(name="x")
+            op.binary_var(name="y")
+            op.binary_var(name="z")
+            op.binary_var(name="w")
+            linear = {"x": 2, "y": 1, "z": -1, "w": 1}
+            quadratic = {("y", "z"): -2, ("w", "w"): 1}
+            op.minimize(linear=linear, quadratic=quadratic)
+
+            linear_constraint = {"x": 1, "w": 1}
+            op.linear_constraint(linear_constraint, Constraint.Sense.LE, 1, "P(xw)")
+            linear_constraint = {"y": 1, "z": 1}
+            op.linear_constraint(linear_constraint, Constraint.Sense.LE, 1, "P(yz)")
+            linear_constraint = {"y": 2, "z": 1}
+            op.linear_constraint(linear_constraint, Constraint.Sense.EQ, 1, "None 1")
+            quadratic_constraint = {("x", "x"): -2, ("y", "w"): 1}
+            op.quadratic_constraint(
+                linear_constraint, quadratic_constraint, Constraint.Sense.LE, 1, "None 2"
+            )
+
+            lip.penalty = 5
+            op2 = lip.convert(op)
+            quadratic[("x", "w")] = lip.penalty
+            quadratic[("y", "z")] = quadratic[("y", "z")] + lip.penalty
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 1)
+            self.assertEqual(op2.get_num_quadratic_constraints(), 1)
+
+    def test_linear_inequality_to_penalty2(self):
+        """Test special constraint to penalty x+y >= 1 -> P(1-x-y+x*y)"""
+
+        op = QuadraticProgram()
+        lip = LinearInequalityToPenalty()
+
+        op.binary_var(name="x")
+        op.binary_var(name="y")
+        # Linear constraints
+        linear_constraint = {"x": 1, "y": 1}
+        op.linear_constraint(linear_constraint, Constraint.Sense.GE, 1, "P(1-x-y+xy)")
+
+        # Test with no max/min
+        with self.subTest("No max/min"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            lip.penalty = 1
+            constant = 1
+            linear = {"x": -1 * lip.penalty, "y": -1 * lip.penalty}
+            quadratic = {("x", "y"): lip.penalty}
+            op2 = lip.convert(op)
+            cnst = op2.objective.constant
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(cnst, constant)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test maximize
+        with self.subTest("Maximize"):
+            linear = {"x": 2, "y": 1}
+            op.maximize(linear=linear)
+            lip.penalty = 5
+            constant = -1 * lip.penalty
+            linear["x"] = linear["x"] + lip.penalty
+            linear["y"] = linear["y"] + lip.penalty
+            quadratic = {("x", "y"): -1 * lip.penalty}
+            op2 = lip.convert(op)
+            cnst = op2.objective.constant
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(cnst, constant)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test minimize
+        with self.subTest("Minimize"):
+            linear = {"x": 2, "y": 1}
+            op.minimize(linear=linear)
+            lip.penalty = 5
+            constant = lip.penalty
+            linear["x"] = linear["x"] - lip.penalty
+            linear["y"] = linear["y"] - lip.penalty
+            quadratic = {("x", "y"): lip.penalty}
+            op2 = lip.convert(op)
+            cnst = op2.objective.constant
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(cnst, constant)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test combination
+        with self.subTest("Combination"):
+            op = QuadraticProgram()
+            lip = LinearInequalityToPenalty()
+
+            op.binary_var(name="x")
+            op.binary_var(name="y")
+            op.binary_var(name="z")
+            op.binary_var(name="w")
+            linear = {"x": 2, "y": 1, "z": -1, "w": 1}
+            quadratic = {("y", "z"): -2, ("w", "w"): 1}
+            op.minimize(linear=linear, quadratic=quadratic)
+
+            linear_constraint = {"x": 1, "w": 1}
+            op.linear_constraint(linear_constraint, Constraint.Sense.GE, 1, "P(1-x-w+xw)")
+            linear_constraint = {"y": 1, "z": 1}
+            op.linear_constraint(linear_constraint, Constraint.Sense.GE, 1, "P(1-y-z+yz)")
+            linear_constraint = {"y": 2, "z": 1}
+            op.linear_constraint(linear_constraint, Constraint.Sense.EQ, 1, "None 1")
+            quadratic_constraint = {("x", "x"): -2, ("y", "w"): 1}
+            op.quadratic_constraint(
+                linear_constraint, quadratic_constraint, Constraint.Sense.LE, 1, "None 2"
+            )
+
+            lip.penalty = 5
+            op2 = lip.convert(op)
+            linear["x"] = linear["x"] - lip.penalty
+            linear["y"] = linear["y"] - lip.penalty
+            linear["z"] = linear["z"] - lip.penalty
+            linear["w"] = linear["w"] - lip.penalty
+            quadratic[("x", "w")] = lip.penalty
+            quadratic[("y", "z")] = quadratic[("y", "z")] + lip.penalty
+            constant = lip.penalty
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(cnst, constant)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 1)
+            self.assertEqual(op2.get_num_quadratic_constraints(), 1)
+
+    def test_linear_inequality_to_penalty3(self):
+        """Test special constraint to penalty x1+x2+x3+... >= n-1 -> P(x1*x2+x1*x3+...)"""
+
+        op = QuadraticProgram()
+        lip = LinearInequalityToPenalty()
+        op.binary_var_list(5)
+
+        # Linear constraints
+        n = 5
+        op.linear_constraint([1, 1, 1, 1, 1], Constraint.Sense.GE, n - 1, "")
+
+        # Test with no max/min
+        with self.subTest("No max/min"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            penalty = 5
+            lip.penalty = penalty
+
+            constant = 10
+            linear = [n - 1, n - 1, n - 1, n - 1, n - 1]
+            quadratic = [
+                [0, 1, 1, 1, 1],
+                [0, 0, 1, 1, 1],
+                [0, 0, 0, 1, 1],
+                [0, 0, 0, 0, 1],
+                [0, 0, 0, 0, 0],
+            ]
+            op2 = lip.convert(op)
+            cnst2 = op2.objective.constant / penalty
+            ldct2 = op2.objective.linear.to_array() / penalty * -1
+            qdct2 = op2.objective.quadratic.to_array() / penalty
+            self.assertEqual(cnst2, constant)
+            self.assertEqual(ldct2.tolist(), linear)
+            self.assertEqual(qdct2.tolist(), quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+    def test_linear_inequality_to_penalty4(self):
+        """Test special constraint to penalty x1+x2+x3+... <= 1 -> P(x1*x2+x1*x3+...)"""
+
+        op = QuadraticProgram()
+        lip = LinearInequalityToPenalty()
+
+        op.binary_var(name="x")
+        op.binary_var(name="y")
+        op.binary_var(name="z")
+        # Linear constraints
+        linear_constraint = {"x": 1, "y": 1, "z": 1}
+        op.linear_constraint(linear_constraint, Constraint.Sense.LE, 1, "P(xy+yz+zx)")
+
+        # Test with no max/min
+        with self.subTest("No max/min"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            penalty = 1
+            quadratic = {("x", "y"): penalty, ("x", "z"): penalty, ("y", "z"): penalty}
+            op2 = lip.convert(op)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test maximize
+        op = QuadraticProgram()
+        op.binary_var_list(5)
+        linear2 = [1, 1, 0, 0, 0]
+        op.maximize(linear=linear2)
+        op.linear_constraint([1, 1, 1, 1, 1], Constraint.Sense.LE, 1, "")
+
+        with self.subTest("Maximum"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            lip.penalty = 5
+            quadratic2 = [
+                [0, 1, 1, 1, 1],
+                [0, 0, 1, 1, 1],
+                [0, 0, 0, 1, 1],
+                [0, 0, 0, 0, 1],
+                [0, 0, 0, 0, 0],
+            ]
+            op2 = lip.convert(op)
+            ldct2 = op2.objective.linear.to_array()
+            qdct2 = op2.objective.quadratic.to_array() / lip.penalty * -1
+            self.assertEqual(ldct2.tolist(), linear2)
+            self.assertEqual(qdct2.tolist(), quadratic2)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+    def test_linear_inequality_to_penalty6(self):
+        """Test special constraint to penalty 6 x-y <= 0 -> P(x-x*y)"""
+
+        op = QuadraticProgram()
+        lip = LinearInequalityToPenalty()
+
+        op.binary_var(name="x")
+        op.binary_var(name="y")
+        # Linear constraints
+        linear_constraint = {"x": 1, "y": -1}
+        op.linear_constraint(linear_constraint, Constraint.Sense.LE, 0, "P(x-xy)")
+
+        # Test with no max/min
+        with self.subTest("No max/min"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            penalty = 1
+            linear = {"x": penalty}
+            quadratic = {("x", "y"): -1 * penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test maximize
+        with self.subTest("Maximize"):
+            linear = {"x": 2, "y": 1}
+            op.maximize(linear=linear)
+            penalty = 4
+            linear["x"] = linear["x"] - penalty
+            quadratic = {("x", "y"): penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test minimize
+        with self.subTest("Minimize"):
+            linear = {"x": 2, "y": 1}
+            op.minimize(linear={"x": 2, "y": 1})
+            penalty = 4
+            linear["x"] = linear["x"] + penalty
+            quadratic = {("x", "y"): -1 * penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+    def test_linear_inequality_to_penalty7(self):
+        """Test special constraint to penalty 6 x-y >= 0 -> P(y-x*y)"""
+
+        op = QuadraticProgram()
+        lip = LinearInequalityToPenalty()
+
+        op.binary_var(name="x")
+        op.binary_var(name="y")
+        # Linear constraints
+        linear_constraint = {"x": 1, "y": -1}
+        op.linear_constraint(linear_constraint, Constraint.Sense.GE, 0, "P(y-xy)")
+
+        # Test with no max/min
+        with self.subTest("No max/min"):
+            self.assertEqual(op.get_num_linear_constraints(), 1)
+            penalty = 1
+            linear = {"y": penalty}
+            quadratic = {("x", "y"): -1 * penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test maximize
+        with self.subTest("Maximize"):
+            linear = {"x": 2, "y": 1}
+            op.maximize(linear=linear)
+            penalty = 4
+            linear["y"] = linear["y"] - penalty
+            quadratic = {("x", "y"): penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+        # Test minimize
+        with self.subTest("Minimize"):
+            linear = {"x": 2, "y": 1}
+            op.minimize(linear={"x": 2, "y": 1})
+            penalty = 4
+            linear["y"] = linear["y"] + penalty
+            quadratic = {("x", "y"): -1 * penalty}
+            op2 = lip.convert(op)
+            ldct = op2.objective.linear.to_dict(use_name=True)
+            qdct = op2.objective.quadratic.to_dict(use_name=True)
+            self.assertEqual(ldct, linear)
+            self.assertEqual(qdct, quadratic)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+    def test_quadratic_program_to_qubo_inequality_to_penalty(self):
+        """Test QuadraticProgramToQubo, passing inequality pattern"""
+
+        op = QuadraticProgram()
+        conv = QuadraticProgramToQubo()
+        op.binary_var(name="x")
+        op.binary_var(name="y")
+
+        # Linear constraints
+        linear_constraint = {"x": 1, "y": 1}
+        op.linear_constraint(linear_constraint, Constraint.Sense.GE, 1, "P(1-x-y+xy)")
+
+        conv.penalty = 1
+        constant = 1
+        linear = {"x": -conv.penalty, "y": -conv.penalty}
+        quadratic = {("x", "y"): conv.penalty}
+        op2 = conv.convert(op)
+        cnst = op2.objective.constant
+        ldct = op2.objective.linear.to_dict(use_name=True)
+        qdct = op2.objective.quadratic.to_dict(use_name=True)
+        self.assertEqual(cnst, constant)
+        self.assertEqual(ldct, linear)
+        self.assertEqual(qdct, quadratic)
+        self.assertEqual(op2.get_num_linear_constraints(), 0)
+
+    def test_inequality_to_penalty_auto_define_penalty(self):
+        """Test _auto_define_penalty() in InequalityToPenalty"""
+        op = QuadraticProgram()
+        op.integer_var(name="x", lowerbound=1, upperbound=3)
+        op.integer_var(name="y", lowerbound=-1, upperbound=4)
+        op.integer_var(name="z", lowerbound=-5, upperbound=-1)
+        op.maximize(linear={"x": 1, "y": 1, "z": 1})
+        lip = LinearInequalityToPenalty()
+        self.assertEqual(lip._auto_define_penalty(op), 12)
+        op = QuadraticProgram()
+        op.integer_var(name="x", lowerbound=1, upperbound=3)
+        op.integer_var(name="y", lowerbound=-1, upperbound=4)
+        op.integer_var(name="z", lowerbound=-5, upperbound=-1)
+        op.maximize(linear={"x": -1, "y": -1, "z": -1})
+        lip = LinearInequalityToPenalty()
+        self.assertEqual(lip._auto_define_penalty(op), 12)
+        op = QuadraticProgram()
+        op.integer_var(name="x", lowerbound=1, upperbound=3)
+        op.integer_var(name="y", lowerbound=-1, upperbound=4)
+        op.integer_var(name="z", lowerbound=-5, upperbound=-1)
+        op.maximize(quadratic={(0, 0): 1, (0, 1): 1, (0, 2): 1, (1, 1): 1, (1, 2): 1, (2, 2): 1})
+        lip = LinearInequalityToPenalty()
+        self.assertEqual(lip._auto_define_penalty(op), 103)
+        op = QuadraticProgram()
+        op.integer_var(name="x", lowerbound=1, upperbound=3)
+        op.integer_var(name="y", lowerbound=-1, upperbound=4)
+        op.integer_var(name="z", lowerbound=-5, upperbound=-1)
+        op.maximize(
+            quadratic={(0, 0): -1, (0, 1): -1, (0, 2): -1, (1, 1): -1, (1, 2): -1, (2, 2): -1}
+        )
+        lip = LinearInequalityToPenalty()
+        self.assertEqual(lip._auto_define_penalty(op), 103)
+        op = QuadraticProgram()
+        op.integer_var(lowerbound=-2, upperbound=1, name="x")
+        op.minimize(quadratic={("x", "x"): 1})
+        lip = LinearInequalityToPenalty()
+        self.assertEqual(lip._auto_define_penalty(op), 5)
+
+    def test_linear_inequality_to_penalty8(self):
+        """Test combinations of inequality constraints"""
+
+        with self.subTest("minimize 1"):
+            op = QuadraticProgram()
+            op.binary_var("x")
+            op.binary_var("y")
+            op.binary_var("z")
+            op.integer_var(-1, 4, "q")
+            op.minimize(linear={"x": 1, "y": 1, "z": 1}, quadratic={("q", "q"): -1})
+            op.linear_constraint({"x": 1, "y": -1}, "<=", 0)
+            op.linear_constraint({"x": 1, "y": 1, "z": 1}, "<=", 1)
+            op2 = LinearInequalityToPenalty().convert(op)
+            self.assertEqual(op2.get_num_vars(), 4)
+            self.assertEqual(op2.get_num_binary_vars(), 3)
+            self.assertEqual(op2.get_num_integer_vars(), 1)
+            self.assertEqual(op2.get_num_continuous_vars(), 0)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+            self.assertEqual(op2.get_num_quadratic_constraints(), 0)
+            obj = op2.objective
+            self.assertEqual(obj.constant, 0)
+            self.assertDictEqual(obj.linear.to_dict(use_name=True), {"x": 21, "y": 1, "z": 1})
+            self.assertDictEqual(
+                obj.quadratic.to_dict(use_name=True),
+                {("x", "z"): 20, ("y", "z"): 20, ("q", "q"): -1},
+            )
+
+        with self.subTest("maximize 1"):
+            op = QuadraticProgram()
+            op.binary_var("x")
+            op.binary_var("y")
+            op.binary_var("z")
+            op.integer_var(-1, 4, "q")
+            op.maximize(linear={"x": 1, "y": 1, "z": 1}, quadratic={("q", "q"): -1})
+            op.linear_constraint({"x": 1, "y": -1}, "<=", 0)
+            op.linear_constraint({"x": 1, "y": 1, "z": 1}, "<=", 1)
+            op2 = LinearInequalityToPenalty().convert(op)
+            self.assertEqual(op2.get_num_vars(), 4)
+            self.assertEqual(op2.get_num_binary_vars(), 3)
+            self.assertEqual(op2.get_num_integer_vars(), 1)
+            self.assertEqual(op2.get_num_continuous_vars(), 0)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+            self.assertEqual(op2.get_num_quadratic_constraints(), 0)
+            obj = op2.objective
+            self.assertEqual(obj.constant, 0)
+            self.assertDictEqual(obj.linear.to_dict(use_name=True), {"x": -19, "y": 1, "z": 1})
+            self.assertDictEqual(
+                obj.quadratic.to_dict(use_name=True),
+                {("x", "z"): -20, ("y", "z"): -20, ("q", "q"): -1},
+            )
+
+        with self.subTest("minimize 2"):
+            op = QuadraticProgram()
+            op.binary_var("x")
+            op.binary_var("y")
+            op.binary_var("z")
+            op.integer_var(-1, 4, "q")
+            op.minimize(linear={"x": 1, "y": 1, "z": 1}, quadratic={("q", "q"): -1})
+            op.linear_constraint({"x": 1, "y": -1}, ">=", 0)
+            op.linear_constraint({"x": 1, "y": 1, "z": 1}, ">=", 2)
+            op2 = LinearInequalityToPenalty().convert(op)
+            self.assertEqual(op2.get_num_vars(), 4)
+            self.assertEqual(op2.get_num_binary_vars(), 3)
+            self.assertEqual(op2.get_num_integer_vars(), 1)
+            self.assertEqual(op2.get_num_continuous_vars(), 0)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+            self.assertEqual(op2.get_num_quadratic_constraints(), 0)
+            obj = op2.objective
+            self.assertEqual(obj.constant, 60)
+            self.assertDictEqual(obj.linear.to_dict(use_name=True), {"x": -39, "y": -19, "z": -39})
+            self.assertDictEqual(
+                obj.quadratic.to_dict(use_name=True),
+                {("x", "z"): 20, ("y", "z"): 20, ("q", "q"): -1},
+            )
+
+        with self.subTest("maximize 2"):
+            op = QuadraticProgram()
+            op.binary_var("x")
+            op.binary_var("y")
+            op.binary_var("z")
+            op.integer_var(-1, 4, "q")
+            op.maximize(linear={"x": 1, "y": 1, "z": 1}, quadratic={("q", "q"): -1})
+            op.linear_constraint({"x": 1, "y": -1}, ">=", 0)
+            op.linear_constraint({"x": 1, "y": 1, "z": 1}, ">=", 2)
+            op2 = LinearInequalityToPenalty().convert(op)
+            self.assertEqual(op2.get_num_vars(), 4)
+            self.assertEqual(op2.get_num_binary_vars(), 3)
+            self.assertEqual(op2.get_num_integer_vars(), 1)
+            self.assertEqual(op2.get_num_continuous_vars(), 0)
+            self.assertEqual(op2.get_num_linear_constraints(), 0)
+            self.assertEqual(op2.get_num_quadratic_constraints(), 0)
+            obj = op2.objective
+            self.assertEqual(obj.constant, -60)
+            self.assertDictEqual(obj.linear.to_dict(use_name=True), {"x": 41, "y": 21, "z": 41})
+            self.assertDictEqual(
+                obj.quadratic.to_dict(use_name=True),
+                {("x", "z"): -20, ("y", "z"): -20, ("q", "q"): -1},
+            )
+
+
+if __name__ == "__main__":
+    unittest.main()