make_symmetric operation

_test/test_poly_matrix.py

@@ -1,5 +1,7 @@
 import numpy as np
 import pytest
+import scipy.sparse as sp
+import scipy.linalg as la
 
 from poly_matrix import PolyMatrix, sorted_dict
 
@@ -448,6 +450,49 @@ def test_multiply():
     np.testing.assert_allclose(S.get_matrix_dense(["x1", "x2"]), S_test)
 
 
+def test_init_from_sparse():
+    """Define polymatrix from sparse matrix. Test symmetrization."""
+    # Get sparse matrix (already symmetric)
+    Q1 = get_Q()
+    var_dict = Q1.generate_variable_dict()
+    assert Q1.symmetric, ValueError("Starting matrix not symmetric. Test broken")
+    # Convert to sparse matrix
+    Q1_sparse = Q1.get_matrix()
+    # Convert back to polymatrix
+    Q2, _ = PolyMatrix.init_from_sparse(Q1_sparse, var_dict)
+    Q2_sparse = Q2.get_matrix()
+    # Test equality
+    np.testing.assert_allclose(
+        Q1_sparse.todense(), Q2_sparse.todense(), err_msg="Matrices not equal"
+    )
+
+    # Make matrix assymmetric by removing one block
+    Q2_sparse[4, 0] = 0.0
+    Q2_sparse[4, 1] = 0.0
+    Q2_sparse.eliminate_zeros()
+    Q2_mat = Q2_sparse.todense()
+    # Define new polymatrices
+    Q3_asym, _ = PolyMatrix.init_from_sparse(Q2_sparse, var_dict)
+    Q3_sym, _ = PolyMatrix.init_from_sparse(Q2_sparse, var_dict, symmetric=True)
+    # Check symmetry
+    Q3_asym_mat = Q3_asym.get_matrix(var_dict).todense()
+    assert not la.issymmetric(Q3_asym_mat), ValueError("Matrix should be assymmetric")
+    Q3_sym_mat = Q3_sym.get_matrix(var_dict).todense()
+    assert la.issymmetric(Q3_sym_mat), ValueError("Matrix should be symmetric")
+    # Check asymmetric matrix mapped properly
+    np.testing.assert_allclose(
+        Q2_mat,
+        Q3_asym_mat,
+        err_msg="Asymmetric matrix definition fail.",
+    )
+    # Check that symmetry and get_matrix operations commute
+    np.testing.assert_allclose(
+        (Q2_mat + Q2_mat.T) / 2,
+        Q3_sym_mat,
+        err_msg="symmetry and matrix retrieval ops don't commute.",
+    )
+
+
 if __name__ == "__main__":
     test_multiply()
 
@@ -465,4 +510,6 @@ def test_multiply():
     test_operations_simple()
     test_operations_advanced()
 
+    test_init_from_sparse()
+
     print("all tests passed")

poly_matrix/poly_matrix.py

@@ -102,8 +102,9 @@ def __init__(self, symmetric=True):
 
         self.symmetric = symmetric
 
-        # dictionary of form {variable-key : variable-size}
+        # dictionary of form {variable-key : variable-size} for enforcing
         # TODO(FD) consider replacing with NamedTuple
+        # consistent variable sizes for matrix blocks.
         self.variable_dict_i = {}
         self.variable_dict_j = {}
 
@@ -126,7 +127,7 @@ def init_from_row_list(row_list, row_labels=None):
         return poly_vstack
 
     @staticmethod
-    def init_from_sparse(A, var_dict, unfold=False):
+    def init_from_sparse(A, var_dict, unfold=False, symmetric=False):
         """Construct polymatrix from sparse matrix (e.g. from learning method)"""
         self = PolyMatrix(symmetric=False)
         var_dict_augmented = augment(var_dict)
@@ -150,10 +151,41 @@ def init_from_sparse(A, var_dict, unfold=False):
         if unfold:
             new_var_dict = {val[2]: 1 for val in var_dict_augmented.values()}
             return self, new_var_dict
+
+        # Symmetrize if required
+        if symmetric:
+            self.make_symmetric()
+
         return self, var_dict
 
-    def drop(self, variables):
-        for v in variables:
+    def make_symmetric(self):
+        """Convert a polymatrix from assymmetric to symmetric.
+        Performed in place."""
+        if self.symmetric:
+            return
+        # Search through keys and identify
+        for key1 in self.matrix.keys():
+            for key2 in self.matrix[key1].keys():
+                # check if flipped order element exists (always holds for diagonals)
+                if key2 in self.matrix.keys() and key1 in self.matrix[key2].keys():
+                    # Symmetrize stored element
+                    mat1 = self.matrix[key1][key2]
+                    mat2 = self.matrix[key2][key1]
+                    mat = (mat1 + mat2.T) / 2
+                    self.matrix[key1][key2] = mat
+                    if not key1 == key2:
+                        self.matrix[key2][key1] = mat.T
+                else:  # If other side not populated, assume its zero
+                    mat = self.matrix[key1][key2] / 2
+                    self.matrix[key1][key2] = mat
+                    self.__setitem__((key2, key1), val=mat.T, symmetric=False)
+        # Align variable list
+        self.variable_dict_j = self.variable_dict_i
+        # Set flag
+        self.symmetric = True
+
+    def drop(self, variables_i):
+        for v in variables_i:
             if v in self.matrix:
                 self.matrix.pop(v)
             if v in self.variable_dict_i:
@@ -640,7 +672,7 @@ def get_block_matrices(self, key_list=None):
                     blocks.append(np.zeros((i_size, j_size)))
         return blocks
 
-    def get_expr(self, variables=None):
+    def get_expr(self, variables=None, homog=None):
 
         if not self.symmetric:
             warnings.warn(
@@ -675,19 +707,27 @@ def get_expr(self, variables=None):
                             if i > 0:
                                 sub_expr += " + "
                             if diag:
-                                sub_expr += (
-                                    "{:.3f}".format(vals[i])
-                                    + "*"
-                                    + ":".join([key1, str(row[i])])
-                                    + "^2"
-                                )
+                                if key1 == homog:
+                                    sub_expr += "{:.3f}".format(vals[i])
+                                else:
+                                    sub_expr += (
+                                        "{:.3f}".format(vals[i])
+                                        + "*"
+                                        + ":".join([key1, str(row[i])])
+                                        + "^2"
+                                    )
                             else:
+                                if key1 == homog:
+                                    key1str = ""
+                                else:
+                                    key1str = "*" + ":".join([key1, str(row[i])])
+                                if key2 == homog:
+                                    key2str = ""
+                                else:
+                                    key2str = "*" + ":".join([key2, str(col[i])])
+
                                 sub_expr += (
-                                    "{:.3f}".format(vals[i] * 2)
-                                    + "*"
-                                    + ":".join([key1, str(row[i])])
-                                    + "*"
-                                    + ":".join([key2, str(col[i])])
+                                    "{:.3f}".format(vals[i] * 2) + key1str + key2str
                                 )
                         if first_expr:
                             expr += sub_expr + "\n"