Merge branch 'master' into use-sparse-solvers

.gitmodules

@@ -1,9 +1,6 @@
 [submodule "poly_matrix"]
 	path = poly_matrix
 	url = [email protected]:utiasASRL/poly_matrix
-[submodule "starloc"]
-	path = starloc
-	url = [email protected]:utiasASRL/starloc
 [submodule "certifiable-tools"]
 	path = certifiable-tools
 	url = [email protected]:utiasASRL/certifiable-tools

examples/robust_lifter.py

@@ -13,13 +13,20 @@ def get_problem(robust=True):
     # Create lifter
     np.random.seed(0)
     n_landmarks = 4
-    d= 3
+    d = 3
     if robust:
         n_outliers = 1
     else:
         n_outliers = 0
 
-    lifter = Lifter(d=d, n_landmarks=n_landmarks+n_outliers, robust=robust, n_outliers=n_outliers, level=level, variable_list=None)
+    lifter = Lifter(
+        d=d,
+        n_landmarks=n_landmarks + n_outliers,
+        robust=robust,
+        n_outliers=n_outliers,
+        level=level,
+        variable_list=None,
+    )
     Q, y = lifter.get_Q()
 
     from auto_template.learner import Learner
@@ -71,4 +78,4 @@ def plot_problem(prob, lifter, fname=""):
         prob, lifter = get_problem(robust=robust)
         fname = f"certifiable-tools/_examples/test_prob_{number}G.pkl"
         save_test_problem(**prob, fname=fname)
-        #plot_problem(prob, lifter, fname=fname.replace(".pkl", ".png"))
+        # plot_problem(prob, lifter, fname=fname.replace(".pkl", ".png"))

lifters/base_class.py

@@ -1,9 +1,14 @@
-from abc import ABC, abstractproperty, abstractmethod
+from abc import ABC, abstractmethod, abstractproperty
 
 import numpy as np
 
 
 class BaseClass(ABC):
+    """
+    This BaseClass is to be used as a template for creating new Lifters.
+    All the listed functionalities need to be defined by inheriting functions.
+    """
+
     LEVELS = ["no"]
     PARAM_LEVELS = ["no"]
     VARIABLE_LIST = ["h"]
@@ -32,71 +37,51 @@ def __init__(
             self.variable_list = self.VARIABLE_LIST
 
         # variables that get overwritten upon initialization
+        self.parameters = [1.0]
         self.theta_ = None
         self.var_dict_ = None
         self.y_ = None
 
         self.d = d
-        self.generate_random_setup()
-
-    @property
-    def d(self):
-        return self.d_
-
-    @d.setter
-    def d(self, var):
-        assert var in [1, 2, 3]
-        self.d_ = var
 
     @abstractproperty
     def var_dict(self):
         """Return key,size pairs of all variables."""
         return
 
-    @abstractmethod
-    def get_param_idx_dict(self, var_subset=None):
-        """Return key,index pairs of all parameters touched by var_subset"""
+    @abstractproperty
+    def theta(self):
+        """Return ground truth theta."""
         return
 
     @abstractmethod
-    def get_level_dims(self, n=1):
+    def get_x(self, theta, var_subset=None) -> np.ndarray:
         return
 
     @abstractmethod
-    def generate_random_setup(self):
+    def get_p(self, parameters=None, var_subset=None) -> np.ndarray:
         return
 
-    # @abstractmethod
-    def generate_random_theta(self):
-        pass
-
-    # @abstractmethod
+    @abstractmethod
     def sample_theta(self):
         return
 
-    # @abstractmethod
-    def sample_parameters(self, x=None):
-        return
-
     @abstractmethod
-    def get_x(self, theta, var_subset=None) -> np.ndarray:
-        return
+    def get_Q(self, noise=1e-3):
+        Warning("get_Q not implemented yet")
+        return None, None
 
-    @abstractmethod
-    def get_p(self, parameters=None, var_subset=None) -> np.ndarray:
-        return
+    def get_A_known(self):
+        return []
+
+    def sample_parameters(self, x=None):
+        if self.param_level == "no":
+            return [1.0]
 
-    # @abstractmethod
     def get_parameters(self, var_subset=None) -> list:
-        return
+        if self.param_level == "no":
+            return [1.0]
 
     def get_grad(self, t, y):
         Warning("get_grad not implement yet")
         return None
-
-    def get_Q(self, noise=1e-3):
-        Warning("get_Q not implemented yet")
-        return None, None
-
-    def get_A_known(self):
-        return []

lifters/state_lifter.py

@@ -48,6 +48,7 @@ def unravel_multi_index_triu(flat_indices, shape):
 
 
 class StateLifter(BaseClass):
+    HOM = "h"
     # consider singular value zero below this
     EPS_SVD = 1e-5
 
@@ -171,10 +172,6 @@ def theta(self):
             self.theta_ = self.generate_random_theta()
         return self.theta_
 
-    @theta.setter
-    def theta(self, value):
-        self.theta_ = value
-
     @property
     def base_var_dict(self):
         var_dict = {"x": self.d**2 + self.d}
@@ -189,7 +186,7 @@ def sub_var_dict(self):
     @property
     def var_dict(self):
         if self.var_dict_ is None:
-            self.var_dict_ = {"h": 1}
+            self.var_dict_ = {self.HOM: 1}
             self.var_dict_.update(self.base_var_dict)
             self.var_dict_.update(self.sub_var_dict)
         return self.var_dict_
@@ -410,7 +407,7 @@ def var_list_row(self, var_subset=None, force_parameters_off=False):
 
         label_list = []
         if force_parameters_off:
-            param_dict = {"h": 0}
+            param_dict = {self.HOM: 0}
         else:
             param_dict = self.get_param_idx_dict(var_subset)
         for idx, key in enumerate(param_dict.keys()):
@@ -451,7 +448,7 @@ def get_basis_from_poly_rows(self, basis_poly_list, var_subset=None):
         for i, bi_poly in enumerate(basis_poly_list):
             # test that this constraint holds
 
-            bi = bi_poly.get_matrix((["h"], all_dict))
+            bi = bi_poly.get_matrix(([self.HOM], all_dict))
 
             if bi.shape[1] == self.get_dim_X(var_subset) * self.get_dim_P():
                 ai = self.get_reduced_a(bi, var_subset=var_subset)
@@ -477,13 +474,15 @@ def get_basis_from_poly_rows(self, basis_poly_list, var_subset=None):
     def get_vector_dense(self, poly_row_sub):
         # complete the missing variables
         var_dict = self.var_dict_row()
-        poly_row_all = poly_row_sub.get_matrix((["h"], var_dict), output_type="poly")
+        poly_row_all = poly_row_sub.get_matrix(
+            ([self.HOM], var_dict), output_type="poly"
+        )
         vector = np.empty(0)
         for param in self.get_param_idx_dict().keys():
             # extract each block corresponding to a bigger matrix
             sub_mat = PolyMatrix(symmetric=True)
             for vari, varj in itertools.combinations_with_replacement(self.var_dict, 2):
-                val = poly_row_all["h", f"{param}.{vari}.{varj}"]
+                val = poly_row_all[self.HOM, f"{param}.{vari}.{varj}"]
                 if np.ndim(val) > 0:
                     if vari != varj:
                         sub_mat[vari, varj] = val.reshape(
@@ -522,13 +521,13 @@ def convert_polyrow_to_Apoly(self, poly_row, correct=True):
             keyi_m, keyj_n = var_keys.split(".")
             m = keyi_m.split(":")[-1]
             n = keyj_n.split(":")[-1]
-            if param in ["h", "h.h"]:
+            if param in [self.HOM, f"{self.HOM}.{self.HOM}"]:
                 param_val = 1.0
             else:
                 param_val = parameters[param_dict[param]]
 
             # divide off-diagonal elements by sqrt(2)
-            newval = poly_row["h", key] * param_val
+            newval = poly_row[self.HOM, key] * param_val
             if correct and not ((keyi_m == keyj_n) and (m == n)):
                 newval /= np.sqrt(2)
 
@@ -581,7 +580,7 @@ def convert_a_to_polyrow(
         for keyi, keyj in itertools.combinations_with_replacement(var_dict, 2):
             if keyi in poly_mat.matrix and keyj in poly_mat.matrix[keyi]:
                 val = poly_mat.matrix[keyi][keyj]
-                labels = self.get_labels("h", keyi, keyj)
+                labels = self.get_labels(self.HOM, keyi, keyj)
                 if keyi != keyj:
                     vals = val.flatten()
                 else:
@@ -590,7 +589,7 @@ def convert_a_to_polyrow(
                 assert len(labels) == len(vals)
                 for label, v in zip(labels, vals):
                     if np.any(np.abs(v) > self.EPS_SPARSE):
-                        poly_row["h", label] = v
+                        poly_row[self.HOM, label] = v
         return poly_row
 
     def convert_b_to_polyrow(self, b, var_subset, tol=1e-10) -> PolyMatrix:
@@ -605,7 +604,7 @@ def convert_b_to_polyrow(self, b, var_subset, tol=1e-10) -> PolyMatrix:
         mask = np.abs(b) > tol
         var_list = [v for i, v in enumerate(self.var_list_row(var_subset)) if mask[i]]
         for key, val in zip(var_list, b[mask]):
-            poly_row["h", key] = val
+            poly_row[self.HOM, key] = val
         return poly_row
 
     def apply_templates(
@@ -727,12 +726,10 @@ def apply_template(
                             raise IndexError(
                                 "something went wrong in augment_basis_list"
                             )
-                except ValueError:
+                except ValueError as e:
                     pass
-                new_poly_row["h", key_ij] = bi_poly["h", key]
+                new_poly_row[self.HOM, key_ij] = bi_poly["h", key]
             new_poly_rows.append(new_poly_row)
-        if verbose:
-            print("done")
         return new_poly_rows
 
     def get_vec_around_gt(self, delta: float = 0):
@@ -788,7 +785,7 @@ def generate_Y(self, factor=FACTOR, ax=None, var_subset=None):
                 else:
                     ax.scatter(*theta[:, :2].T)
 
-            x = self.get_x(theta, parameters, var_subset=var_subset)
+            x = self.get_x(theta=theta, parameters=parameters, var_subset=var_subset)
             X = np.outer(x, x)
 
             # generates [1*x, a1*x, ..., aK*x]
@@ -857,7 +854,7 @@ def get_reduced_a(self, bi, var_subset=None, sparse=False):
         if isinstance(bi, np.ndarray):
             len_b = len(bi)
         elif isinstance(bi, PolyMatrix):
-            bi = bi.get_matrix((["h"], self.var_dict_row(var_subset)))
+            bi = bi.get_matrix(([self.HOM], self.var_dict_row(var_subset)))
             len_b = bi.shape[1]
         else:
             # bi can be a scipy sparse matrix,
@@ -996,7 +993,7 @@ def test_constraints(self, A_list, errors: str = "raise", n_seeds: int = 3):
                 np.random.seed(i)
                 t = self.sample_theta()
                 p = self.get_parameters()
-                x = self.get_x(t, p)
+                x = self.get_x(theta=t, parameters=p)
 
                 constraint_violation = abs(x.T @ A @ x)
                 max_violation = max(max_violation, constraint_violation)
@@ -1019,7 +1016,7 @@ def get_A0(self, var_subset=None):
         else:
             var_dict = self.var_dict
         A0 = PolyMatrix()
-        A0["h", "h"] = 1.0
+        A0[self.HOM, self.HOM] = 1.0
         return A0.get_matrix(var_dict)
 
     def get_A_b_list(self, A_list, var_subset=None):
@@ -1037,20 +1034,22 @@ def get_param_idx_dict(self, var_subset=None):
         - if param_level == 'ppT': {'l': 0, 'p_0:0.p_0:0': 1, ..., 'p_0:d-1:.p_0:d-1': 1}
         """
         if self.param_level == "no":
-            return {"h": 0}
+            return {self.HOM: 0}
 
         if var_subset is None:
             var_subset = self.var_dict
         variables = self.get_variable_indices(var_subset)
-        param_keys = ["h"] + [f"p_{i}:{d}" for i in variables for d in range(self.d)]
+        param_keys = [self.HOM] + [
+            f"p_{i}:{d}" for i in variables for d in range(self.d)
+        ]
         if self.param_level == "p":
             param_dict = {p: i for i, p in enumerate(param_keys)}
         elif self.param_level == "ppT":
             i = 0
             param_dict = {}
             for pi, pj in itertools.combinations_with_replacement(param_keys, 2):
-                if pi == pj == "h":
-                    param_dict["h"] = i
+                if pi == pj == self.HOM:
+                    param_dict[self.HOM] = i
                 else:
                     param_dict[f"{pi}.{pj}"] = i
                 i += 1

lifters/stereo_lifter.py

@@ -448,7 +448,7 @@ def get_error(self, xtheta_hat):
 
     def local_solver_manopt(self, t0, y, W=None, verbose=False, method="CG", **kwargs):
         import pymanopt
-        from pymanopt.manifolds import SpecialOrthogonalGroup, Euclidean, Product
+        from pymanopt.manifolds import Euclidean, Product, SpecialOrthogonalGroup
 
         if method == "CG":
             from pymanopt.optimizers import ConjugateGradient as Optimizer  # fastest