Uniform penalty method (#227)

Roughly, the `Instance.penalty_method` introduced in #200 converts a
constrained problem

$$
\begin{align*}
  \min_x &f(x) & \\
  \text{s.t.} &g_i(x) = 0 & (\forall i)
\end{align*}
$$

into an unconstrained parametric problem with parameters $\lambda_i$

$$
\min_x f(x) + \sum_i \lambda_i g_i(x)^2
$$

On this PR, `Instance.uniform_penalty_method` is introduced. This
converts above into a single parameter unconstrained problem

$$
\min_x f(x) + \lambda \sum_i g_i(x)^2
$$

python/ommx/ommx/_ommx_rust.pyi

@@ -105,7 +105,9 @@ class Instance:
     def validate(self) -> None: ...
     def as_pubo_format(self) -> dict: ...
     def as_qubo_format(self) -> tuple[dict, float]: ...
+    def as_parametric_instance(self) -> ParametricInstance: ...
     def penalty_method(self) -> ParametricInstance: ...
+    def uniform_penalty_method(self) -> ParametricInstance: ...
 
 class Linear:
     @staticmethod

python/ommx/ommx/v1/__init__.py

@@ -421,12 +421,176 @@ def as_pubo_format(self) -> dict[tuple[int, ...], float]:
         return instance.as_pubo_format()
 
     def penalty_method(self) -> ParametricInstance:
-        """
-        Convert the instance to a parametric instance for penalty method.
+        r"""
+        Convert to a parametric unconstrained instance by penalty method.
+
+        Roughly, this converts a constrained problem
+
+        .. math::
+
+            \begin{align*}
+                \min_x & \space f(x) & \\
+                \text{ s.t. } & \space g_i(x) = 0 & (\forall i)
+            \end{align*}
+
+        to an unconstrained problem with parameters
+
+        .. math::
+
+            \min_x f(x) + \sum_i \lambda_i g_i(x)^2
+
+        where :math:`\lambda_i` are the penalty weight parameters for each constraint.
+        If you want to use single weight parameter, use :py:meth:`uniform_penalty_method` instead.
+
+        The removed constrains are stored in :py:attr:`~ParametricInstance.removed_constraints`.
+
+        :raises RuntimeError: If the instance contains inequality constraints.
+
+        Examples
+        =========
+
+        .. doctest::
+
+            >>> from ommx.v1 import Instance, DecisionVariable, Constraint
+            >>> x = [DecisionVariable.binary(i) for i in range(3)]
+            >>> instance = Instance.from_components(
+            ...     decision_variables=x,
+            ...     objective=sum(x),
+            ...     constraints=[x[0] + x[1] == 1, x[1] + x[2] == 1],
+            ...     sense=Instance.MAXIMIZE,
+            ... )
+            >>> instance.objective
+            Function(x0 + x1 + x2)
+
+            >>> pi = instance.penalty_method()
+
+            The constraint is put in `removed_constraints`
+
+            >>> pi.get_constraints()
+            []
+            >>> len(pi.get_removed_constraints())
+            2
+            >>> pi.get_removed_constraints()[0]
+            RemovedConstraint(Function(x0 + x1 - 1) == 0, reason=penalty_method)
+            >>> pi.get_removed_constraints()[1]
+            RemovedConstraint(Function(x1 + x2 - 1) == 0, reason=penalty_method)
+
+            There are two parameters corresponding to the two constraints
+
+            >>> len(pi.get_parameters())
+            2
+            >>> p1 = pi.get_parameters()[0]
+            >>> p1.id, p1.name
+            (3, 'penalty_weight')
+            >>> p2 = pi.get_parameters()[1]
+            >>> p2.id, p2.name
+            (4, 'penalty_weight')
+
+            Substitute all parameters to zero to get the original objective
+
+            >>> instance0 = pi.with_parameters({p1.id: 0.0, p2.id: 0.0})
+            >>> instance0.objective
+            Function(x0 + x1 + x2)
+
+            Substitute all parameters to one
+
+            >>> instance1 = pi.with_parameters({p1.id: 1.0, p2.id: 1.0})
+            >>> instance1.objective
+            Function(x0*x0 + 2*x0*x1 + 2*x1*x1 + 2*x1*x2 + x2*x2 - x0 - 3*x1 - x2 + 2)
+
         """
         instance = _ommx_rust.Instance.from_bytes(self.to_bytes())
         return ParametricInstance.from_bytes(instance.penalty_method().to_bytes())
 
+    def uniform_penalty_method(self) -> ParametricInstance:
+        r"""
+        Convert to a parametric unconstrained instance by penalty method with uniform weight.
+
+        Roughly, this converts a constrained problem
+
+        .. math::
+
+            \begin{align*}
+                \min_x & \space f(x) & \\
+                \text{ s.t. } & \space g_i(x) = 0 & (\forall i)
+            \end{align*}
+
+        to an unconstrained problem with a parameter
+
+        .. math::
+
+            \min_x f(x) + \lambda \sum_i g_i(x)^2
+
+        where :math:`\lambda` is the uniform penalty weight parameter for all constraints.
+
+        The removed constrains are stored in :py:attr:`~ParametricInstance.removed_constraints`.
+
+        :raises RuntimeError: If the instance contains inequality constraints.
+
+        Examples
+        =========
+
+        .. doctest::
+
+            >>> from ommx.v1 import Instance, DecisionVariable
+            >>> x = [DecisionVariable.binary(i) for i in range(3)]
+            >>> instance = Instance.from_components(
+            ...     decision_variables=x,
+            ...     objective=sum(x),
+            ...     constraints=[sum(x) == 3],
+            ...     sense=Instance.MAXIMIZE,
+            ... )
+            >>> instance.objective
+            Function(x0 + x1 + x2)
+
+            >>> pi = instance.uniform_penalty_method()
+
+            The constraint is put in `removed_constraints`
+
+            >>> pi.get_constraints()
+            []
+            >>> len(pi.get_removed_constraints())
+            1
+            >>> pi.get_removed_constraints()[0]
+            RemovedConstraint(Function(x0 + x1 + x2 - 3) == 0, reason=uniform_penalty_method)
+
+            There is only one parameter in the instance
+
+            >>> len(pi.get_parameters())
+            1
+            >>> p = pi.get_parameters()[0]
+            >>> p.id
+            3
+            >>> p.name
+            'uniform_penalty_weight'
+
+            Substitute `p = 0` to get the original objective
+
+            >>> instance0 = pi.with_parameters({p.id: 0.0})
+            >>> instance0.objective
+            Function(x0 + x1 + x2)
+
+            Substitute `p = 1`
+
+            >>> instance1 = pi.with_parameters({p.id: 1.0})
+            >>> instance1.objective
+            Function(x0*x0 + 2*x0*x1 + 2*x0*x2 + x1*x1 + 2*x1*x2 + x2*x2 - 5*x0 - 5*x1 - 5*x2 + 9)
+
+        """
+        instance = _ommx_rust.Instance.from_bytes(self.to_bytes())
+        return ParametricInstance.from_bytes(
+            instance.uniform_penalty_method().to_bytes()
+        )
+
+    def as_parametric_instance(self) -> ParametricInstance:
+        """
+        Convert the instance to a :class:`ParametricInstance`.
+        """
+        instance = _ommx_rust.Instance.from_bytes(self.to_bytes())
+        return ParametricInstance.from_bytes(
+            instance.as_parametric_instance().to_bytes()
+        )
+
 
 @dataclass
 class ParametricInstance(InstanceBase):
@@ -2028,10 +2192,31 @@ class RemovedConstraint:
 
     raw: _RemovedConstraint
 
+    def __repr__(self) -> str:
+        reason = f"reason={self.removed_reason}"
+        if self.removed_reason_parameters:
+            reason += ", " + ", ".join(
+                f"{key}={value}"
+                for key, value in self.removed_reason_parameters.items()
+            )
+        if self.equality == Equality.EQUALITY_EQUAL_TO_ZERO:
+            return f"RemovedConstraint({self.function.__repr__()} == 0, {reason})"
+        if self.equality == Equality.EQUALITY_LESS_THAN_OR_EQUAL_TO_ZERO:
+            return f"RemovedConstraint({self.function.__repr__()} <= 0, {reason})"
+        return self.raw.__repr__()
+
+    @property
+    def equality(self) -> Equality.ValueType:
+        return self.raw.constraint.equality
+
     @property
     def id(self) -> int:
         return self.raw.constraint.id
 
+    @property
+    def function(self) -> Function:
+        return Function(self.raw.constraint.function)
+
     @property
     def name(self) -> str | None:
         return (

python/ommx/src/instance.rs

@@ -37,8 +37,16 @@ impl Instance {
         Ok((serde_pyobject::to_pyobject(py, &qubo)?.extract()?, constant))
     }
 
-    pub fn penalty_method(&self) -> ParametricInstance {
-        ParametricInstance(self.0.clone().penalty_method())
+    pub fn as_parametric_instance(&self) -> ParametricInstance {
+        ParametricInstance(self.0.clone().into())
+    }
+
+    pub fn penalty_method(&self) -> Result<ParametricInstance> {
+        Ok(ParametricInstance(self.0.clone().penalty_method()?))
+    }
+
+    pub fn uniform_penalty_method(&self) -> Result<ParametricInstance> {
+        Ok(ParametricInstance(self.0.clone().uniform_penalty_method()?))
     }
 }
 

rust/ommx/src/convert/instance.rs

@@ -1,7 +1,7 @@
 use crate::v1::{
     decision_variable::Kind,
     instance::{Description, Sense},
-    Function, Instance, Parameter, ParametricInstance, RemovedConstraint,
+    Equality, Function, Instance, Parameter, ParametricInstance, RemovedConstraint,
 };
 use anyhow::{bail, Context, Result};
 use approx::AbsDiffEq;
@@ -135,15 +135,18 @@ impl Instance {
             .boxed()
     }
 
-    pub fn penalty_method(self) -> ParametricInstance {
+    pub fn penalty_method(self) -> Result<ParametricInstance> {
         let id_base = self.defined_ids().last().map(|id| id + 1).unwrap_or(0);
         let mut objective = self.objective().into_owned();
         let mut parameters = Vec::new();
         let mut removed_constraints = Vec::new();
         for (i, c) in self.constraints.into_iter().enumerate() {
+            if c.equality() != Equality::EqualToZero {
+                bail!("Penalty method is only for equality constraints. Non-equality constraint is found: ID={}", c.id);
+            }
             let parameter = Parameter {
                 id: id_base + i as u64,
-                name: Some("penalty".to_string()),
+                name: Some("penalty_weight".to_string()),
                 subscripts: vec![c.id as i64],
                 ..Default::default()
             };
@@ -156,7 +159,7 @@ impl Instance {
                 removed_reason_parameters: Default::default(),
             });
         }
-        ParametricInstance {
+        Ok(ParametricInstance {
             description: self.description,
             objective: Some(objective),
             constraints: Vec::new(),
@@ -165,7 +168,42 @@ impl Instance {
             parameters,
             constraint_hints: self.constraint_hints,
             removed_constraints,
+        })
+    }
+
+    pub fn uniform_penalty_method(self) -> Result<ParametricInstance> {
+        let id_base = self.defined_ids().last().map(|id| id + 1).unwrap_or(0);
+        let mut objective = self.objective().into_owned();
+        let parameter = Parameter {
+            id: id_base,
+            name: Some("uniform_penalty_weight".to_string()),
+            ..Default::default()
+        };
+        let mut removed_constraints = Vec::new();
+        let mut quad_sum = Function::zero();
+        for c in self.constraints.into_iter() {
+            if c.equality() != Equality::EqualToZero {
+                bail!("Uniform penalty method is only for equality constraints. Non-equality constraint is found: ID={}", c.id);
+            }
+            let f = c.function().into_owned();
+            quad_sum = quad_sum + f.clone() * f;
+            removed_constraints.push(RemovedConstraint {
+                constraint: Some(c),
+                removed_reason: "uniform_penalty_method".to_string(),
+                removed_reason_parameters: Default::default(),
+            });
         }
+        objective = objective + &parameter * quad_sum;
+        Ok(ParametricInstance {
+            description: self.description,
+            objective: Some(objective),
+            constraints: Vec::new(),
+            decision_variables: self.decision_variables.clone(),
+            sense: self.sense,
+            parameters: vec![parameter],
+            constraint_hints: self.constraint_hints,
+            removed_constraints,
+        })
     }
 
     pub fn binary_ids(&self) -> BTreeSet<u64> {
@@ -474,10 +512,43 @@ mod tests {
 
         #[test]
         fn test_penalty_method(instance in Instance::arbitrary()) {
-            let parametric_instance = instance.clone().penalty_method();
+            let Ok(parametric_instance) = instance.clone().penalty_method() else { return Ok(()); };
+            let dv_ids = parametric_instance.defined_decision_variable_ids();
+            let p_ids = parametric_instance.defined_parameter_ids();
+            prop_assert!(dv_ids.is_disjoint(&p_ids));
+
+            let used_ids = parametric_instance.used_ids().unwrap();
+            let all_ids = dv_ids.union(&p_ids).cloned().collect();
+            prop_assert!(used_ids.is_subset(&all_ids));
+
+            // Put every penalty weights to zero
+            let parameters = Parameters {
+                entries: p_ids.iter().map(|&id| (id, 0.0)).collect(),
+            };
+            let substituted = parametric_instance.clone().with_parameters(parameters).unwrap();
+            prop_assert!(instance.objective().abs_diff_eq(&substituted.objective(), 1e-10));
+            prop_assert_eq!(substituted.constraints.len(), 0);
+
+            // Put every penalty weights to two
+            let parameters = Parameters {
+                entries: p_ids.iter().map(|&id| (id, 2.0)).collect(),
+            };
+            let substituted = parametric_instance.with_parameters(parameters).unwrap();
+            let mut objective = instance.objective().into_owned();
+            for c in &instance.constraints {
+                let f = c.function().into_owned();
+                objective = objective + 2.0 * f.clone() * f;
+            }
+            prop_assert!(objective.abs_diff_eq(&substituted.objective(), 1e-10));
+        }
+
+        #[test]
+        fn test_uniform_penalty_method(instance in Instance::arbitrary()) {
+            let Ok(parametric_instance) = instance.clone().uniform_penalty_method() else { return Ok(()); };
             let dv_ids = parametric_instance.defined_decision_variable_ids();
             let p_ids = parametric_instance.defined_parameter_ids();
             prop_assert!(dv_ids.is_disjoint(&p_ids));
+            prop_assert_eq!(p_ids.len(), 1);
 
             let used_ids = parametric_instance.used_ids().unwrap();
             let all_ids = dv_ids.union(&p_ids).cloned().collect();