ldlt inverse now verified for matrices up to and including size 3x3

.gitignore

@@ -116,6 +116,7 @@ TODO.txt
 b_asic/_version.py
 docs_sphinx/_build/
 docs_sphinx/examples
+docs_sphinx/sg_execution_times.rst
 result_images/
 .coverage
 Digraph.gv

b_asic/core_operations.py

@@ -1123,7 +1123,7 @@ class MADS(AbstractOperation):
 
     def __init__(
         self,
-        is_add: Optional[bool] = False,
+        is_add: Optional[bool] = True,
         src0: Optional[SignalSourceProvider] = None,
         src1: Optional[SignalSourceProvider] = None,
         src2: Optional[SignalSourceProvider] = None,
@@ -1142,15 +1142,24 @@ class MADS(AbstractOperation):
             latency_offsets=latency_offsets,
             execution_time=execution_time,
         )
-        self._is_add = is_add
-        # self.set_param("is_add", is_add)
+        self.set_param("is_add", is_add)
 
     @classmethod
     def type_name(cls) -> TypeName:
         return TypeName("mads")
 
     def evaluate(self, a, b, c):
-        return a + b * c if self._is_add else a - b * c
+        return a + b * c if self.is_add else a - b * c
+
+    @property
+    def is_add(self) -> bool:
+        """Get if operation is an addition."""
+        return self.param("is_add")
+
+    @is_add.setter
+    def is_add(self, is_add: bool) -> None:
+        """Set if operation is an addition."""
+        self.set_param("is_add", is_add)
 
     @property
     def is_linear(self) -> bool:

b_asic/sfg_generators.py

@@ -456,40 +456,41 @@ def ldlt_matrix_inverse(N: int, is_complex: bool) -> SFG:
 
     # R*di*R^T factorization
     for i in range(N):
-        for k in range(1, i):
-            D[i] = MADS(False, D[i], M[k][i], R[k][i])
+        for k in range(i):
+            D[i] = MADS(False, D[i], M[k][i], R[k][i], name="M1")
 
         D_inv[i] = Reciprocal(D[i])
 
-        for j in range(i, N):
+        for j in range(i + 1, N):
             R[i][j] = A[i][j]
 
-            for k in range(1, i):
-                R[i][j] = MADS(False, R[i][j], M[k][i], R[k][j])
+            for k in range(i):
+                R[i][j] = MADS(False, R[i][j], M[k][i], R[k][j], name="M2")
 
             if is_complex:
                 M[i][j] = ComplexConjugate(R[i][j])
             else:
                 M[i][j] = R[i][j]
 
-            R[i][j] = MADS(True, Constant(0, name="0"), R[i][j], D_inv[i])
+            R[i][j] = MADS(True, Constant(0, name="0"), R[i][j], D_inv[i], name="M3")
 
     # back substitution
     A_inv = [[None for _ in range(N)] for _ in range(N)]
     for i in reversed(range(N)):
         A_inv[i][i] = D_inv[i]
-        for j in reversed(range(i)):
-            for k in reversed(range(j, N)):
+        for j in reversed(range(i + 1)):
+            for k in reversed(range(j + 1, N)):
                 if k == N - 1 and i != j:
-                    # my custom
-                    if A_inv[k][i]:
+                    A_inv[j][i] = MADS(
+                        False, Constant(0, name="0"), R[j][k], A_inv[i][k], name="M4"
+                    )
+                else:
+                    if A_inv[i][k]:
                         A_inv[j][i] = MADS(
-                            False, Constant(0, name="0"), R[j][k], A_inv[k][i]
+                            False, A_inv[j][i], R[j][k], A_inv[i][k], name="M5"
                         )
                     else:
-                        A_inv[j][i] = Constant(0, name="0")
-                else:
-                    A_inv[j][i] = MADS(False, A_inv[j][i], R[j][k], A_inv[k][i])
+                        A_inv[j][i] = MADS(False, A_inv[j][i], R[j][k], A_inv[k][i])
 
     outputs = []
     for i in range(N):

test/test_sfg_generators.py

@@ -12,6 +12,7 @@ from b_asic.sfg_generators import (
     direct_form_1_iir,
     direct_form_2_iir,
     direct_form_fir,
+    ldlt_matrix_inverse,
     radix_2_dif_fft,
     transposed_direct_form_fir,
     wdf_allpass,
@@ -637,3 +638,58 @@ class TestRadix2FFT:
         POINTS = 5
         with pytest.raises(ValueError, match="Points must be a power of two."):
             radix_2_dif_fft(points=POINTS)
+
+
+class TestLdltMatrixInverse:
+    def test_1x1(self):
+        sfg = ldlt_matrix_inverse(N=1, is_complex=False)
+
+        A_input = [Constant(5)]
+
+        sim = Simulation(sfg, A_input)
+        sim.run_for(1)
+
+        res = sim.results
+        assert np.isclose(res["0"], 0.2)
+
+    def test_2x2_simple_spd(self):
+        sfg = ldlt_matrix_inverse(N=2, is_complex=False)
+
+        A = np.array([[1, 2], [2, 1]])
+        A_input = [Constant(1), Constant(2), Constant(1)]
+
+        A_inv = np.linalg.inv(A)
+
+        sim = Simulation(sfg, A_input)
+        sim.run_for(1)
+
+        res = sim.results
+        assert np.isclose(res["0"], A_inv[0, 0])
+        assert np.isclose(res["1"], A_inv[0, 1])
+        assert np.isclose(res["2"], A_inv[1, 1])
+
+    def test_3x3_simple_spd(self):
+        sfg = ldlt_matrix_inverse(N=3, is_complex=False)
+
+        A = np.array([[2, -1, 0], [-1, 3, -1], [0, -1, 2]])
+        A_input = [
+            Constant(2),
+            Constant(-1),
+            Constant(0),
+            Constant(3),
+            Constant(-1),
+            Constant(2),
+        ]
+
+        A_inv = np.linalg.inv(A)
+
+        sim = Simulation(sfg, A_input)
+        sim.run_for(1)
+
+        res = sim.results
+        assert np.isclose(res["0"], A_inv[0, 0])
+        assert np.isclose(res["1"], A_inv[0, 1])
+        assert np.isclose(res["2"], A_inv[0, 2])
+        assert np.isclose(res["3"], A_inv[1, 1])
+        assert np.isclose(res["4"], A_inv[1, 2])
+        assert np.isclose(res["5"], A_inv[2, 2])