Add certificate de-symmetrization

This has made me realize that lexorder is different in InflationLP and InflationSDP. We may want to correct this

Author: apozas

inflation/lp/InflationLP.py

@@ -873,6 +873,49 @@ def evaluate_certificate(self, prob_array: np.ndarray) -> float:
                  "valid for other distributions.")
         return self.evaluate_polynomial(self.certificate_as_dict(), prob_array)
 
+    def desymmetrize_certificate(self) -> dict:
+        """If the scenario contains symmetries other than the inflation
+        symmetries, this function writes a certificate of infeasibility valid
+        for non-symmetric distributions too.
+
+        Parameters
+        ----------
+        replace : bool, optional
+            Whether to replace the certificate in ``self.solution_object``.
+            By default ``False``.
+
+        Returns
+        -------
+        dict
+            The expression of the un-symmetrized certificate in terms of
+            probabilities and marginals. The certificate of incompatibility is
+            ``cert < 0``.
+        """
+        try:
+            dual = self.solution_object["dual_certificate"]
+        except AttributeError:
+            raise Exception("For extracting a certificate you need to solve " +
+                            "a problem. Call \"InflationSDP.solve()\" first.")
+
+        desymmetrized = {}
+        norm = len(self.InflationProblem.symmetries)
+        lexmon_names = self.InflationProblem._lexrepr_to_copy_index_free_names
+        for symm in self.InflationProblem.symmetries:
+            for mon, coeff in dual.items():
+                mon = self.monomial_from_name[mon]
+                if not mon.is_zero:
+                    desymm_mon = lexmon_names[symm[mon.as_lexmon]]
+                    desymm_mon = sorted(desymm_mon, key=lambda x: x[0])
+                    if not mon.is_one:
+                        desymm_name = "P[" + " ".join(desymm_mon) + "]"
+                    else:
+                        desymm_name = "1"
+                    if desymm_name not in desymmetrized:
+                        desymmetrized[desymm_name] = coeff / norm
+                    else:
+                        desymmetrized[desymm_name] += coeff / norm
+        return desymmetrized
+
     ###########################################################################
     # OTHER ROUTINES EXPOSED TO THE USER                                      #
     ##########################################################################

inflation/sdp/InflationSDP.py

@@ -1121,6 +1121,43 @@ def evaluate_certificate(self, prob_array: np.ndarray) -> float:
                  "valid for other distributions.")
         return self.evaluate_polynomial(self.certificate_as_dict(), prob_array)
 
+    def desymmetrize_certificate(self) -> dict:
+        """If the scenario contains symmetries other than the inflation
+        symmetries, this function writes a certificate of infeasibility valid
+        for non-symmetric distributions too.
+
+        Returns
+        -------
+        dict
+            The expression of the un-symmetrized certificate in terms of
+            probabilities and marginals. The certificate of incompatibility is
+            ``cert < 0``.
+        """
+        try:
+            dual = self.solution_object["dual_certificate"]
+        except AttributeError:
+            raise Exception("For extracting a certificate you need to solve " +
+                            "a problem. Call \"InflationSDP.solve()\" first.")
+
+        desymmetrized = {}
+        norm = len(self.InflationProblem.symmetries)
+        lexmon_names = self.InflationProblem._lexrepr_to_copy_index_free_names
+        for symm in self.InflationProblem.symmetries:
+            for mon, coeff in dual.items():
+                mon = self.monomial_from_name[mon]
+                if not mon.is_zero:
+                    desymm_mon = lexmon_names[symm[mon.as_lexmon-1]]
+                    desymm_mon = sorted(desymm_mon, key=lambda x: x[0])
+                    if not mon.is_one:
+                        desymm_name = "P[" + " ".join(desymm_mon) + "]"
+                    else:
+                        desymm_name = "1"
+                    if desymm_name not in desymmetrized:
+                        desymmetrized[desymm_name] = coeff / norm
+                    else:
+                        desymmetrized[desymm_name] += coeff / norm
+        return desymmetrized
+
     ###########################################################################
     # OTHER ROUTINES EXPOSED TO THE USER                                      #
     ###########################################################################

test/test_symmetry.py

@@ -2,7 +2,7 @@
 import warnings
 import numpy as np
 
-from inflation import InflationProblem
+from inflation import InflationProblem, InflationLP
 from inflation.symmetry_utils import (discover_distribution_symmetries,
                                       group_elements_from_generators)
 
@@ -24,17 +24,19 @@ def setUpClass(cls):
     bellScenario = InflationProblem({"Lambda": ["A", "B"]},
                                     outcomes_per_party=[2, 2],
                                     settings_per_party=[2, 2],
-                                    inflation_level_per_source=[1])
+                                    inflation_level_per_source=[1],
+                                    classical_sources='all')
 
     triangle = InflationProblem({"Lambda": ["A", "B"],
                                  "Mu": ["B", "C"],
                                  "Sigma": ["C", "A"]},
                                 outcomes_per_party=[2, 2, 2],
                                 inflation_level_per_source=[2, 1, 1])
 
+    PRbox_symmetries = discover_distribution_symmetries(PR_box,
+                                                        bellScenario)
+
     def test_discover(self):
-        PRbox_symmetries = discover_distribution_symmetries(self.PR_box,
-                                                            self.bellScenario)
         # Order: (a=0,x=0), (a=1,x=0), (a=0,x=1), (a=1,x=1),
         #        (b=0,y=0), (b=1,y=0), (b=0,y=1), (b=1,y=1)
         # There's five symmetries: flip parties, flip X, flip Y, and flip the
@@ -60,7 +62,7 @@ def test_discover(self):
         swapped = [symm[4:] + symm[:4] for symm in symmetries]
         symmetries = symmetries + swapped
         self.assertSetEqual(set(map(tuple, symmetries)),
-                            set(map(tuple, PRbox_symmetries)),
+                            set(map(tuple, self.PRbox_symmetries)),
                             "Failed to discover the symmetries of the PR box.")
 
     def test_discover_inflation(self):
@@ -89,3 +91,24 @@ def test_group_elements_from_generators(self):
         self.assertSetEqual(set(map(tuple, elements)),
                             set(map(tuple, truth)),
                             "Failed to generate S3 from generators.")
+
+    def test_desymmetrized_certificate(self):
+        self.bellScenario.add_symmetries(self.PRbox_symmetries)
+        lp = InflationLP(self.bellScenario, verbose=0)
+        lp.set_distribution(self.PR_box)
+        lp.solve()
+        certificate = lp.desymmetrize_certificate()
+        truth = {
+            'P[A_0=0]': 0.125, 'P[A_0=1]': 0.125, 'P[A_1=0]': 0.125, 'P[A_1=1]': 0.125,
+            'P[B_0=0]': 0.125, 'P[B_0=1]': 0.125, 'P[B_1=0]': 0.125, 'P[B_1=1]': 0.125,
+            'P[A_0=0 B_0=0]': -0.1875, 'P[A_0=0 B_0=1]': 0.0625,
+            'P[A_0=0 B_1=0]': -0.1875, 'P[A_0=0 B_1=1]': 0.0625,
+            'P[A_0=1 B_1=1]': -0.1875, 'P[A_0=1 B_1=0]': 0.0625,
+            'P[A_0=1 B_0=1]': -0.1875, 'P[A_0=1 B_0=0]': 0.0625,
+            'P[A_1=0 B_0=0]': -0.1875, 'P[A_1=0 B_0=1]': 0.0625,
+            'P[A_1=0 B_1=1]': -0.1875, 'P[A_1=0 B_1=0]': 0.0625,
+            'P[A_1=1 B_1=0]': -0.1875, 'P[A_1=1 B_1=1]': 0.0625,
+            'P[A_1=1 B_0=1]': -0.1875, 'P[A_1=1 B_0=0]': 0.0625}
+        self.assertDictEqual(certificate, truth,
+                             "Failed to desymmetrize the CHSH inequality.")
+