Fix Homotopy generator using Groebner basis

python/examples/repetition_code.py

@@ -69,51 +69,59 @@ from hccd import path_tracker, homotopy_generator
 
 
 def track_path(args):
-    H = np.array([[1, 1, 1]])
+    H = np.array([[1, 1, 0, 0],
+                  [0, 1, 1, 0],
+                  [0, 0, 1, 1]])
 
-    # TODO: Make this work
     homotopy = homotopy_generator.HomotopyGenerator(H)
 
+    print(f"G: {homotopy.G}")
+    print(f"F: {homotopy.F}")
+    print(f"H: {homotopy.H}")
+    print(f"DH: {homotopy.DH}")
+
     tracker = path_tracker.PathTracker(homotopy, args.euler_step_size, args.euler_max_tries,
                                        args.newton_max_iter, args.newton_convergence_threshold, args.sigma)
 
     ys_start, ys_prime, ys_hat_e, ys = [], [], [], []
 
-    try:
+    y = np.zeros(5) + np.array([0.5, 0.48, -0.1, 0.2, 0 ])
-        y = np.zeros(3)
+    for i in range(args.num_iterations):
-        for i in range(args.num_iterations):
+        y_start, y_prime, y_hat_e, y = tracker.transparent_step(y)
-            y_start, y_prime, y_hat_e, y = tracker.transparent_step(y)
+        ys_start.append(y_start)
-            ys_start.append(y_start)
+        ys_prime.append(y_prime)
-            ys_prime.append(y_prime)
+        ys_hat_e.append(y_hat_e)
-            ys_hat_e.append(y_hat_e)
+        ys.append(y)
-            ys.append(y)
+        print(f"Iteration {i}: {y}")
-            print(f"Iteration {i}: {y}")
-    except Exception as e:
-        print(f"Error: {e}")
 
     ys_start = np.array(ys_start)
     ys_prime = np.array(ys_prime)
     ys_hat_e = np.array(ys_hat_e)
     ys = np.array(ys)
 
-    df = pd.DataFrame({"x1b": ys_start[:, 0],
+    for y in np.transpose(ys):
-                       "x2b": ys_start[:, 1],
+        plt.plot(y)
-                       "tb":  ys_start[:, 2],
-                       "x1p": ys_prime[:, 0],
-                       "x2p": ys_prime[:, 1],
-                       "tp":  ys_prime[:, 2],
-                       "x1e": ys_hat_e[:, 0],
-                       "x2e": ys_hat_e[:, 1],
-                       "te":  ys_hat_e[:, 2],
-                       "x1n": ys[:, 0],
-                       "x2n": ys[:, 1],
-                       "tn":  ys[:, 2]
-                       })
 
-    if args.output:
+    plt.show()
-        df.to_csv(args.output, index=False)
+
-    else:
+    # df = pd.DataFrame({"x1b": ys_start[:, 0],
-        print(df)
+    #                    "x2b": ys_start[:, 1],
+    #                    "tb":  ys_start[:, 2],
+    #                    "x1p": ys_prime[:, 0],
+    #                    "x2p": ys_prime[:, 1],
+    #                    "tp":  ys_prime[:, 2],
+    #                    "x1e": ys_hat_e[:, 0],
+    #                    "x2e": ys_hat_e[:, 1],
+    #                    "te":  ys_hat_e[:, 2],
+    #                    "x1n": ys[:, 0],
+    #                    "x2n": ys[:, 1],
+    #                    "tn":  ys[:, 2]
+    #                    })
+    #
+    # if args.output:
+    #     df.to_csv(args.output, index=False)
+    # else:
+    #     print(df)
 
 
 def main():

python/hccd/homotopy_generator.py

@@ -17,78 +17,63 @@ class HomotopyGenerator:
         self.H_matrix = parity_check_matrix
         self.num_checks, self.num_vars = parity_check_matrix.shape
 
-        # Create symbolic variables
         self.x_vars = [sp.symbols(f'x{i+1}') for i in range(self.num_vars)]
         self.t = sp.symbols('t')
 
-        # Generate G, F, and H
         self.G = self._create_G()
         self.F = self._create_F()
         self.H = self._create_H()
+        self.DH = self._create_DH(self.H)
 
-        # Convert to callable functions
         self._H_lambda = self._create_H_lambda()
         self._DH_lambda = self._create_DH_lambda()
 
     def _create_G(self) -> List[sp.Expr]:
-        """Create G polynomial system (the starting system)."""
-        # For each variable xi, add the polynomial [xi]
         G = []
         for var in self.x_vars:
             G.append(var)
 
         return G
 
-    def _create_F(self) -> List[sp.Expr]:
+    def _create_F(self) -> sp.MutableMatrix:
-        """Create F polynomial system (the target system)."""
         F = []
 
-        # Add 1 - xi^2 for each variable
         for var in self.x_vars:
             F.append(1 - var**2)
 
-        # Add parity check polynomials: 1 - x1*x2*...*xk for each parity check
         for row in self.H_matrix:
-            # Create product of variables that participate in this check
             term = 1
             for i, bit in enumerate(row):
                 if bit == 1:
                     term *= self.x_vars[i]
 
-            if term != 1:  # Only add if there are variables in this check
+            F.append(1 - term)
-                F.append(1 - term)
 
-        return F
+        groebner_basis = sp.groebner(F)
+        return sp.MutableMatrix(groebner_basis)
 
     def _create_H(self) -> List[sp.Expr]:
-        """Create the homotopy H = (1-t)*G + t*F."""
         H = []
 
-        # Make sure G and F have the same length
+        for g, f in zip(self.G, self.F):
-        # Repeat variables from G if needed to match F's length
-        G_extended = self.G.copy()
-        while len(G_extended) < len(self.F):
-            # Cycle through variables to repeat
-            for i in range(min(self.num_vars, len(self.F) - len(G_extended))):
-                G_extended.append(self.x_vars[i % self.num_vars])
-
-        # Create the homotopy
-        for g, f in zip(G_extended, self.F):
             H.append((1 - self.t) * g + self.t * f)
 
         return H
 
+    def _create_DH(self, H: List[sp.Expr]) -> sp.MutableMatrix:
+        all_vars = self.x_vars + [self.t]
+        DH = sp.Matrix([[sp.diff(expr, var)
+                         for var in all_vars] for expr in self.H])
+
+        return DH
+
     def _create_H_lambda(self) -> Callable:
-        """Create a lambda function to evaluate H."""
         all_vars = self.x_vars + [self.t]
         return sp.lambdify(all_vars, self.H, 'numpy')
 
     def _create_DH_lambda(self) -> Callable:
-        """Create a lambda function to evaluate the Jacobian of H."""
         all_vars = self.x_vars + [self.t]
-        jacobian = sp.Matrix([[sp.diff(expr, var)
+        return sp.lambdify(all_vars, self.DH, 'numpy')
-                             for var in all_vars] for expr in self.H])
-        return sp.lambdify(all_vars, jacobian, 'numpy')
 
     def evaluate_H(self, y: np.ndarray) -> np.ndarray:
         """