Merge pull request #164 from ecboghiu/pedro_symmetries

Finish merging the Pedro_Symmetries branch

Author: eliewolfe

inflation/InflationProblem.py

@@ -543,7 +543,7 @@ def _lexrepr_to_names(self) -> np.ndarray:
         return np.asarray([self._interpretation_to_name(
             op_dict,
             include_copy_indices=self._any_inflation)
-                for op_dict in self._lexrepr_to_dicts.flat])
+                for op_dict in self._lexrepr_to_dicts.flat], dtype=object)
 
     @cached_property
     def _original_event_names(self) -> np.ndarray:
@@ -558,7 +558,7 @@ def _original_event_names(self) -> np.ndarray:
         return np.asarray([self._interpretation_to_name(
             self._interpret_operator(event),
             include_copy_indices=False)
-                for event in self.original_dag_events])
+                for event in self.original_dag_events], dtype=object)
 
     @cached_property
     def _lexrepr_to_copy_index_free_names(self) -> np.ndarray:
@@ -576,7 +576,7 @@ def _lexrepr_to_copy_index_free_names(self) -> np.ndarray:
             return np.asarray([self._interpretation_to_name(
                 op_dict,
                 include_copy_indices=False)
-                for op_dict in self._lexrepr_to_dicts.flat])
+                for op_dict in self._lexrepr_to_dicts.flat], dtype=object)
 
     @cached_property
     def _lexrepr_to_all_names(self) -> np.ndarray:
@@ -599,7 +599,7 @@ def _lexrepr_to_all_names(self) -> np.ndarray:
             self._lexrepr_to_copy_index_free_names,
             old_names_v1,
             old_names_v2
-        ), 1)
+        ), 1).astype(object)
 
     @cached_property
     def _lexrepr_to_symbols(self) -> np.ndarray:

inflation/lp/InflationLP.py

@@ -169,13 +169,13 @@ def __init__(self,
             self._boolvec_for_FR_eqs = self.blank_bool_vec
 
         if self.verbose > 1:
-            print("Number of single operator measurements per party:", end="")
+            eprint("Number of single operator measurements per party:", end="")
             prefix = " "
             for i, measures in enumerate(inflationproblem.measurements):
                 op_count = np.prod(measures.shape[:2])
-                print(prefix + f"{self.names[i]}={op_count}", end="")
+                eprint(prefix + f"{self.names[i]}={op_count}", end="")
                 prefix = ", "
-            print()
+            eprint()
         self.use_lpi_constraints = False
 
         self.identity_operator = np.empty((0, self._nr_properties),
@@ -972,7 +972,7 @@ def write_to_file(self, filename: str) -> None:
         # Write file according to the extension
         args = self._prepare_solver_arguments(separate_bounds=True)
         if self.verbose > 0:
-            print("Writing the LP program to", filename)
+            eprint("Writing the LP program to", filename)
         if extension == "lp":
             write_to_lp(args, filename)
         elif extension == "mps":
@@ -1142,7 +1142,7 @@ def _sanitise_monomial(self, mon: Any) -> CompoundMoment:
             try:
                 return self.monomial_from_name[mon]
             except KeyError:
-                print(f"As of now we only recognize \n{list(self.monomial_from_name.keys())}")
+                eprint(f"As of now we only recognize \n{list(self.monomial_from_name.keys())}")
                 return self._sanitise_monomial(self._interpret_name(mon))
         elif isinstance(mon, Real):
             if np.isclose(float(mon), 1):
@@ -1321,7 +1321,7 @@ def _generate_lp(self) -> None:
                 orbits_non_CG, return_index=True, return_inverse=True)
             self.num_non_CG = len(old_reps_non_CG)
             if self.verbose > 1:
-                print(f"Orbits discovered! {self.num_CG} unique monomials.")
+                eprint(f"Orbits discovered! {self.num_CG} unique monomials.")
             # Obtain the real generating monomials after accounting for symmetry
         else:
             self.num_CG = self.raw_n_columns
@@ -1406,7 +1406,7 @@ def _generate_lp(self) -> None:
 
         self._lp_has_been_generated = True
         if self.verbose > 1:
-            print("LP initialization complete, ready to accept further specifics.")
+            eprint("LP initialization complete, ready to accept further specifics.")
 
     def _template_to_event_boolarray(self, template: List[int], decompressor: List[np.ndarray]) -> np.ndarray:
         if template:

inflation/sdp/InflationSDP.py

@@ -35,7 +35,7 @@
                            write_to_mat,
                            write_to_sdpa)
 from ..lp.numbafied import nb_outer_bitwise_or
-from ..utils import clean_coefficients, partsextractor
+from ..utils import clean_coefficients, partsextractor, eprint
 
 
 class InflationSDP:
@@ -101,17 +101,17 @@ def __init__(self,
 
         self.measurements = self._generate_parties()
         if self.verbose > 1:
-            print("Number of single operator measurements per party:", end="")
+            eprint("Number of single operator measurements per party:", end="")
             prefix = " "
             for i, measures in enumerate(self.measurements):
                 counter = count()
                 deque(zip(chain.from_iterable(
                     chain.from_iterable(measures)),
                           counter),
                       maxlen=0)
-                print(prefix + f"{self.names[i]}={next(counter)}", end="")
+                eprint(prefix + f"{self.names[i]}={next(counter)}", end="")
                 prefix = ", "
-            print()
+            eprint()
         self.use_lpi_constraints = False
         self.network_scenario    = inflationproblem.is_network
         self._is_knowable_q_non_networks = \
@@ -149,13 +149,30 @@ def __init__(self,
         self._lexorder = self._default_lexorder.copy()
         self.op_to_lexrepr_dict = {tuple(op): i for i, op in enumerate(self._lexorder)}
         self._lexorder_len = len(self._lexorder)
-        self.lexorder_symmetries = \
+        self.raw_lexorder_symmetries = \
             np.pad(inflationproblem.symmetries + 1, ((0, 0), (1, 0)))
+        # self.lexorder_symmetries = self.raw_lexorder_symmetries.copy()
+        CG_ops = []
+        for boolarray in self._CG_limited_ortho_groups_as_boolarrays:
+            CG_ops.extend(np.flatnonzero(boolarray)+1)
+        CG_ops = np.sort(CG_ops).astype(int)
+        self.lexorder_symmetries=np.array([
+            perm for perm in self.raw_lexorder_symmetries
+            if np.array_equal(np.sort(perm[CG_ops]), CG_ops)
+        ], dtype=int)
+        if self.verbose > 0:
+            old_group_size = len(self.raw_lexorder_symmetries)
+            new_group_size = len(self.lexorder_symmetries)
+            if new_group_size < old_group_size:
+                eprint("Warning: The use of Collins-Gisin notation internally via the argument `include_all_outcomes=False`")
+                eprint(f" means that not all symmetries of the problem can be exploited. Group size drop from {old_group_size} to {new_group_size}.")
 
         self._lexrepr_to_names = \
-            np.hstack((["0"], inflationproblem._lexrepr_to_names))
+            np.hstack((["0"], inflationproblem._lexrepr_to_names)).astype(object)
+        # eprint("CG stuff:", self._lexrepr_to_names[CG_ops])
+        # eprint("else: ", np.setdiff1d(self._lexrepr_to_names, self._lexrepr_to_names[CG_ops]))
         self._lexrepr_to_copy_index_free_names = \
-            np.hstack((["0"], inflationproblem._lexrepr_to_copy_index_free_names))
+            np.hstack((["0"], inflationproblem._lexrepr_to_copy_index_free_names)).astype(object)
         self.op_from_name = {"0": 0}
         for i, op_names in enumerate(inflationproblem._lexrepr_to_all_names.tolist()):
             for op_name in op_names:
@@ -311,7 +328,7 @@ def generate_relaxation(self,
         self.build_columns(column_specification, **kwargs)
         collect()
         if self.verbose > 0:
-            print("Number of columns in the moment matrix:", self.n_columns)
+            eprint("Number of columns in the moment matrix:", self.n_columns)
 
         # Calculate the moment matrix without the inflation symmetries
         unsymmetrized_mm, unsymmetrized_corresp = \
@@ -323,7 +340,7 @@ def generate_relaxation(self,
                          else "")
             if 0 in unsymmetrized_mm.flat:
                 additional_var = 1
-            print("Number of variables" + extra_msg + ":",
+            eprint("Number of variables" + extra_msg + ":",
                   len(unsymmetrized_corresp) + additional_var)
 
         # Calculate the inflation symmetries
@@ -340,7 +357,7 @@ def generate_relaxation(self,
         if self.verbose > 0:
             extra_msg = (" after symmetrization" if symmetrization_required
                          else "")
-            print(f"Number of variables{extra_msg}: "
+            eprint(f"Number of variables{extra_msg}: "
                   + f"{len(self.symmetrized_corresp)+additional_var}")
         del unsymmetrized_mm, unsymmetrized_corresp, \
             symmetrization_required, additional_var
@@ -392,7 +409,7 @@ def generate_relaxation(self,
         self.first_free_idx = first_free_index
         if self.n_vars < old_num_vars:
             if self.verbose > 0:
-                print("Further variable reduction has been made possible. Number of variables in the SDP:",
+                eprint("Further variable reduction has been made possible. Number of variables in the SDP:",
                        self.n_vars)
         # self.compmoment_from_idx = dict(zip(range(self.n_vars), monomials_as_list))
         # self.compmoment_to_idx = dict(zip(monomials_as_list, range(self.n_vars)))
@@ -418,7 +435,7 @@ def generate_relaxation(self,
         self.n_something_knowable = _counter["Semi"]
         self.n_unknowable         = _counter["Unknowable"]
         if self.verbose > 1:
-            print(f"The problem has {self.n_knowable} knowable moments, " +
+            eprint(f"The problem has {self.n_knowable} knowable moments, " +
                   f"{self.n_something_knowable} semi-knowable moments, " +
                   f"and {self.n_unknowable} unknowable moments.")
 
@@ -428,7 +445,7 @@ def generate_relaxation(self,
             self.hermitian_moments = [mon for mon in self.moments
                                       if mon.is_hermitian]
             if self.verbose > 1:
-                print(f"The problem has {len(self.hermitian_moments)} " +
+                eprint(f"The problem has {len(self.hermitian_moments)} " +
                       "non-negative moments.")
 
         # This dictionary useful for certificates_as_probs
@@ -446,7 +463,7 @@ def generate_relaxation(self,
                                                 self.momentmatrix,
                                                 self.verbose)
             if self.verbose > 1 and len(self.idx_level_equalities):
-                print("Number of normalization equalities:",
+                eprint("Number of normalization equalities:",
                       len(self.idx_level_equalities))
             for (norm_idx, summation_idxs) in self.idx_level_equalities:
                 eq_dict = {self.compmoment_from_idx[norm_idx]: 1}
@@ -1404,7 +1421,7 @@ def write_to_file(self, filename: str) -> None:
 
         # Write file according to the extension
         if self.verbose > 0:
-            print("Writing the SDP program to", filename)
+            eprint("Writing the SDP program to", filename)
         if extension == "dat-s":
             write_to_sdpa(self, filename)
         elif extension == "csv":
@@ -1857,7 +1874,7 @@ def _build_cols_from_specs(self, col_specs: List[List[int]]) -> List:
             for specs in col_specs:
                 to_print.append("1" if specs == []
                                 else "".join([self.names[p] for p in specs]))
-            print("Column structure:", "+".join(to_print))
+            eprint("Column structure:", "+".join(to_print))
 
         _zero_lexorder = np.array([0], dtype=np.intc)
         columns      = []
@@ -1969,19 +1986,23 @@ def _discover_columns_symmetries(self) -> np.ndarray:
         permutation_failed = False
         for inf_sym in self.lexorder_symmetries[1:]:
             skip_this_one = False
-            try:
-                total_perm = np.empty(self.n_columns, dtype=int)
-                for i, lexmon in enumerate(self.generating_monomials_1d):
-                    new_lexmon = inf_sym[lexmon]
-                    new_lexmon_canon = self._to_canonical_memoized_1d(
-                        new_lexmon,
-                        apply_only_commutations=True)
+            total_perm = np.empty(self.n_columns, dtype=int)
+            for i, lexmon in enumerate(self.generating_monomials_1d):
+                new_lexmon = np.argsort(inf_sym)[lexmon]
+                new_lexmon_canon = self._to_canonical_memoized_1d(
+                    new_lexmon,
+                    apply_only_commutations=True)
+                try:
                     total_perm[i] \
-                        = self.genmon_1d_to_index[tuple(new_lexmon_canon)]
-            except KeyError:
-                permutation_failed = True
-                permutations_failed += 1
-                skip_this_one = True
+                    = self.genmon_1d_to_index[tuple(new_lexmon_canon)]
+                except KeyError:
+                    eprint(f"Warning: generating monomial before symmetry becomes unrecognizable after symmetry!")
+                    eprint(f" Generating monomial before symmetry: {self._lexrepr_to_names[lexmon]}")
+                    eprint(f" Generating monomial after symmetry: {self._lexrepr_to_names[new_lexmon_canon]}")
+                    permutation_failed = True
+                    permutations_failed += 1
+                    skip_this_one = True
+                    break
             if not skip_this_one:
                 discovered_symmetries.append(total_perm)
         if permutation_failed and (self.verbose > 0):
@@ -2071,14 +2092,14 @@ def _cleanup_after_set_values(self) -> None:
         if self.momentmatrix_has_a_one:
             num_nontrivial_known -= 1
         if self.verbose > 1 and num_nontrivial_known > 0:
-            print("Number of variables with fixed numeric value:",
+            eprint("Number of variables with fixed numeric value:",
                   len(self.known_moments))
         if len(self.semiknown_moments):
             for k in self.known_moments.keys():
                 self.semiknown_moments.pop(k, None)
         num_semiknown = len(self.semiknown_moments)
         if self.verbose > 1 and num_semiknown > 0:
-            print(f"Number of semiknown variables: {num_semiknown}")
+            eprint(f"Number of semiknown variables: {num_semiknown}")
 
     def _reset_lowerbounds(self) -> None:
         """Reset the list of lower bounds."""