Make the simulator, matrices and vectors keep k and stop mutating (#84, #82, #87)

cli/simulator.py

@@ -22,7 +22,7 @@
 parser.add_argument('-c', '--complex_representation', type=str, help='complex representation to use', required=False, choices=["FiveTuple", "Classic"], default="FiveTuple")
 parser.add_argument('-T', '--targets', type=int, help='number of target qubits', required=False, default=1)
 parser.add_argument('-A', '--ancillas', type=int, help='number of ancilla qubits', required=False, default=1)
-parser.add_argument('-m', '--meas', type=int, help='determine measurement outcome - if not set, default is 0 (|0>)', required=False, default=0)
+parser.add_argument('-m', '--meas', type=int, help='determine measurement outcome - if not set, default is 0 (|0>)', required=False, choices=[0, 1], default=0)
 args = parser.parse_args()
 
 if args.complex_representation == "FiveTuple":

src/complex/classic.py

@@ -87,7 +87,8 @@ def __mul__(self, other):
     def __div__(self, other):
         if self.gen is None:
             if isinstance(other, Complex):
-                pass # never used
+                c2_d2 = other.real**2 + other.imag**2
+                return Complex(a=(self.real * other.real + self.imag * other.imag) / c2_d2, b=(self.imag * other.real - self.real * other.imag) / c2_d2)
             elif isinstance(other, float) or isinstance(other, int):
                 return Complex(a=self.real / other, b=self.imag / other, generator=self.gen)
             else:

src/complex/fivetuple.py

@@ -101,7 +101,7 @@ def __mul__(self, other):
             )
 
     def __div__(self, other):
-        if isinstance(other, float):
+        if isinstance(other, float) or isinstance(other, int):
             if self.gen is None:
                 return FiveTuple(a=self.a / other, b=self.b / other, c=self.c / other, d=self.d / other)
             else:

src/complex/matrix.py

@@ -42,7 +42,7 @@ def __init__(self, matrix: list = None, q=None, name=None, generator=None, eleme
                             self.matrix[i][j] = element_representation.one(generator)
                         else:
                             self.matrix[i][j] = element_representation.zero(generator)
-                self.k = 0
+                self.k = k
                 return
 
             self.gen = generator
@@ -100,11 +100,18 @@ def __mul__(self, other):
             return new_matrix
 
         if isinstance(other, float) or isinstance(other, int) or isinstance(other, Complex) or isinstance(other, FiveTuple):
-            # pointwise multiplication
+            # pointwise multiplication into a new matrix
+            new_matrix = self.copy()
             for i in range(len(self.matrix)):
                 for j in range(len(self.matrix[i])):
-                    self.matrix[i][j] = self.matrix[i][j] * other
-            return self
+                    if self.element_representation == FiveTuple and not isinstance(other, FiveTuple):
+                        # five-tuple coefficients are integers, other scales go through k
+                        if not float(other).is_integer():
+                            raise ValueError("a five-tuple matrix can be multiplied only by an integer")
+                        new_matrix.matrix[i][j] = self.matrix[i][j].multiply_by_real(int(other))
+                    else:
+                        new_matrix.matrix[i][j] = self.matrix[i][j] * other
+            return new_matrix
         else:
             raise ValueError("matrix multiplication not supported for this representation")
 
@@ -114,7 +121,7 @@ def __div__(self, other : float) -> "Matrix":
         return self * (1 / other)
 
     def __truediv__(self, other : float) -> "Matrix":
-        return self * (1 / other)
+        return self.__div__(other)
 
     def __repr__(self):
         matrix_str = ""
@@ -123,10 +130,16 @@ def __repr__(self):
             for j in range(len(self.matrix[i])):
                 matrix_str += f"{self.matrix[i][j]}, "
             matrix_str += "]\n"
-        return matrix_str
+        return matrix_str + f"k={self.k}"
 
+    def _check_same_k(self, other : "Matrix"):
+        # five-tuple elements are scaled by 1/sqrt(2)^k, so they can be added only with the same k
+        if self.element_representation == FiveTuple and self.gen is None and self.k != other.k:
+            raise ValueError(f"five-tuple matrices with different k ({self.k} and {other.k}) cannot be added")
+
     def __add__(self, other : "Matrix") -> "Matrix":
         if isinstance(other, Matrix):
+            self._check_same_k(other)
             new_matrix = Matrix(q=len(self.matrix), generator=self.gen, element_representation=self.element_representation, k=self.k)
             for i in range(len(self.matrix)):
                 for j in range(len(other.matrix)):
@@ -137,6 +150,7 @@ def __add__(self, other : "Matrix") -> "Matrix":
 
     def __sub__(self, other : "Matrix") -> "Matrix":
         if isinstance(other, Matrix):
+            self._check_same_k(other)
             new_matrix = Matrix(q=len(self.matrix), generator=self.gen, element_representation=self.element_representation, k=self.k)
             for i in range(len(self.matrix)):
                 for j in range(len(other.matrix)):
@@ -157,10 +171,11 @@ def transpose(self) -> "Matrix":
         return new_matrix
 
     def conjugate(self) -> "Matrix":
+        new_matrix = self.copy()
         for i in range(len(self.matrix)):
             for j in range(len(self.matrix[i])):
-                self.matrix[i][j] = self.matrix[i][j].conjugate()
-        return self
+                new_matrix.matrix[i][j] = self.matrix[i][j].conjugate()
+        return new_matrix
 
     def tensor(self, other : "Matrix") -> "Matrix":
         # U \otimes V
@@ -397,7 +412,8 @@ def density_from_vector(cls, vector: Vector) -> "Matrix":
         for i in range(n):
             for j in range(n):
                 matrix[i][j] = vector[i] * vector[j].conjugate()
-        return cls(matrix=matrix)
+        # FiveTuple: each product of two elements scaled by 1/sqrt(2)^k is scaled by 1/sqrt(2)^(2k)
+        return cls(matrix=matrix, k=2 * vector.k if vector.element_representation == FiveTuple else 0)
 
     def max_value(self) -> int:
         max = 0

src/complex/vector.py

@@ -5,6 +5,7 @@
 @brief: implementation of vector abstraction over the custom complex number representations
 """
 
+import copy
 import numpy as np
 from .classic import Complex
 from .fivetuple import FiveTuple
@@ -53,7 +54,9 @@ def __setitem__(self, i, value):
 
     def __eq__(self, other):
         if self.gen is None:
-            return all(self.vec[i] == other.vec[i] for i in range(len(self.vec))) and self.k == other.k
+            # value equality: a five-tuple vector with k equals its sqrt(2)-scaled copy with k + 1
+            a, b = self.rescale_with(other) if self.element_representation == FiveTuple else (self, other)
+            return all(a.vec[i] == b.vec[i] for i in range(len(a.vec)))
         else:
             eqs = [self.vec[i] == other.vec[i] for i in range(len(self.vec))]
             if self.k is not None and other.k is not None:
@@ -64,10 +67,9 @@ def __repr__(self):
         return f"[{', '.join(str(v) for v in self.vec)}], k={self.k}"
 
     def copy(self):
-        new_vec = Vector(q=len(self.vec), generator=self.gen, element_representation=self.element_representation, k=self.k)
-        for i in range(len(self.vec)):
-            new_vec[i] = self.vec[i].copy()
-        new_vec.k = self.k
+        # keeps name, bound, n and k, copies only the elements
+        new_vec = copy.copy(self)
+        new_vec.vec = [elem.copy() for elem in self.vec]
         return new_vec
 
     def to_real(self, max_k=None):
@@ -182,6 +184,8 @@ def norm(self) -> float:
 
     def measure(self, q : list[int], result : int) -> "Vector":
         # measure qubit q to result r in {0, 1}
+        if result not in (0, 1):
+            raise ValueError(f"a measurement outcome can be only 0 or 1, not {result}")
         if len(q) == 0:
             return self
         new_vec = self.copy()

src/sim.py

@@ -32,30 +32,38 @@ def __init__(self, qasm_file: str = None, matrix: np.array = None, complex_repre
             raise ValueError("matrix or qasm_file has to be provided")
         if matrix is not None and complex_representation != Complex:
             raise ValueError("matrix can be provided only for Classic complex representation")
+        if meas not in (0, 1):
+            raise ValueError("the measurement outcome (meas) can be only 0 or 1")
         self.qasm_file = qasm_file
         self.complex_representation = complex_representation
         self.matrix = matrix
+        self.measurement_outcome = meas
+        self.basis = basis
+        self.no_measurement = no_measurement
+        self.reset_stats()
+
+    def reset_stats(self) -> None:
+        # every simulate_* parses the file again, and parsing starts from here, so the stats describe one simulation
         self.stats = {}
         self.stats['gate_set'] = GateSet()
         self.stats['d'] = 0
         self.stats['q'] = 0
         self.stats['qreg'] = ''
         self.stats['max_k'] = 0
         self.stats['input_circuit'] = []
+        # each measured qubit once, in the order of its first measurement (empty with no_measurement)
         self.stats['measured_qubits'] = []
         # end-of-circuit measurements as (qubit, classical bit or None) and the classical declarations, both in the
         # input's order and OpenQASM 3 syntax, to write them to the synthesized circuit as they are in the input
         self.stats['measurements'] = []
         self.stats['classical_declarations'] = []
-        self.measurement_outcome = meas
-        self.basis = basis
-        self.no_measurement = no_measurement
 
     """
     parse the input qasm file into a sequence of gates with a respective list of qubits
     also collects information about the input circuits -- qubits, gate set ...
     """
     def parse_file(self) -> tuple[list[tuple[str, list[int]]], list[Vector]]:
+        self.reset_stats()
         vectors = []
         if self.qasm_file is None:
             # matrix 2**q x 2**q
@@ -125,6 +133,7 @@ def parse_file(self) -> tuple[list[tuple[str, list[int]]], list[Vector]]:
                     gates.append(('measure', [i]))
                     self.stats['input_circuit'].append(('measure', [i]))
                     self.stats['measurements'].append((i, None if register is None else f"{register}[{i}]"))
+                    self.add_measured_qubits([i])
                 continue
             # find gates using the register name
             qreg_pattern = re.compile(rf'{re.escape(f"{qreg_name}")}\s*\[')
@@ -148,6 +157,8 @@ def parse_file(self) -> tuple[list[tuple[str, list[int]]], list[Vector]]:
                 self.stats['input_circuit'].append((gate, qubits))
                 if gate == 'measure':
                     self.stats['measurements'].append((qubits[0], self.measurement_target(line)))
+                if gate in ('measure', 'meas'):
+                    self.add_measured_qubits(qubits)
                 self.stats['d'] += 1 if gate != 'measure' and gate != 'meas' else 0
                 if gate not in self.stats['gate_set'] and (gate != 'measure' and gate != 'meas'):
                     self.stats['gate_set'].add_gate(gate, qubits=len(qubits))
@@ -156,6 +167,14 @@ def parse_file(self) -> tuple[list[tuple[str, list[int]]], list[Vector]]:
             raise Exception(f"input file parsing failed")
         return gates, vectors
 
+    def add_measured_qubits(self, qubits: list[int]) -> None:
+        # simulate() ignores measurements with no_measurement, so they measure no qubit
+        if self.no_measurement:
+            return
+        for q in qubits:
+            if q not in self.stats['measured_qubits']:
+                self.stats['measured_qubits'].append(q)
+
     def measurement_target(self, line: str) -> str:
         # the classical bit (or register) of `c = measure q`, `bit c = measure q` or `measure q -> c`, None for a plain
         # `measure q`; a declaration in the statement is moved to the classical declarations
@@ -341,7 +360,6 @@ def simulate(self, vectors: list[Vector], parsed_file: list[tuple[str, list[int]
                                 new_vec = vector.copy()
                                 continue
                             new_vec = vector.measure(qubits, self.measurement_outcome)
-                            self.stats['measured_qubits'].extend(qubits)
                             break
                         else:
                             raise NotImplementedError(f"Gate {op} is not yet implemented")
@@ -474,21 +492,12 @@ def simulate_jamiolkowski(self) -> list[tuple[Vector, Vector]]:
     simulate the given input circuit by creating the unitary matrix instead of vectors
     """
     def simulate_matrix(self) -> Matrix:
-        gates, _ = self.parse_file()
-        matrix = Matrix.i(element_representation=self.complex_representation, q=0, qubits=self.stats['q'])
-        for (op, qubits) in gates:
-            if len(qubits) == 1:
-                mat = getattr(Matrix, op)(element_representation=self.complex_representation, q=qubits[0], qubits=self.stats['q'])
-                matrix = mat * matrix
-            elif len(qubits) == 2:
-                mat = getattr(Matrix, op)(element_representation=self.complex_representation, q1=qubits[0], q2=qubits[1], qubits=self.stats['q'])
-                matrix = mat * matrix
-            elif len(qubits) == 3:
-                mat = getattr(Matrix, op)(element_representation=self.complex_representation, q1=qubits[0], q2=qubits[1], q3=qubits[2], qubits=self.stats['q'])
-                matrix = mat * matrix
-            else:
-                raise NotImplementedError(f"Gate {op} with {len(qubits)} qubits is not yet implemented")
-        return matrix
+        gates, vectors = self.parse_file()
+        if not self.no_measurement and any(gate == 'measure' for gate, _ in gates):
+            raise ValueError("a circuit with measurements has no unitary matrix (no_measurement ignores them)")
+        outputs = [output for _, output in self.simulate(vectors, gates)]
+        # every column went through the same gates, so they share k
+        return Matrix(matrix=[[column[row] for column in outputs] for row in range(len(outputs))], k=outputs[0].k)
 
     """
     return the statistics about the input circuit

tests/test_matrix.py

@@ -89,6 +89,35 @@ def test_scalar_operations():
     np.testing.assert_allclose(to_numpy(Matrix.x().__div__(4)), X / 4)
 
 
+def test_scalar_operations_return_a_new_matrix():
+    # #82, F60: H * 2 and S.conjugate() used to overwrite their operand
+    h, s = Matrix.h(), Matrix.s()
+    assert (h * 2) is not h and s.conjugate() is not s
+    np.testing.assert_allclose(to_numpy(h), H, atol=1e-9)
+    np.testing.assert_allclose(to_numpy(s), np.diag([1, 1j]), atol=1e-9)
+
+
+def test_five_tuple_scalar_operations():
+    # #82, F60: FiveTuple * int used to raise AttributeError
+    h = Matrix.h(FiveTuple)
+    np.testing.assert_allclose(to_numpy(h * 2), 2 * H, atol=1e-9)
+    np.testing.assert_allclose(to_numpy(h), H, atol=1e-9)
+    with pytest.raises(ValueError, match="only by an integer"):
+        h / 2
+
+
+def test_five_tuple_matrices_keep_k():
+    # #82, F26: transpose, dagger, + and - of a matrix with k = 2 used to return k = 0
+    u = Matrix.h(FiveTuple) * Matrix.t(FiveTuple) * Matrix.h(FiveTuple)
+    expected = H @ np.diag([1, OMEGA]) @ H
+    np.testing.assert_allclose(to_numpy(u.transpose()), expected.T, atol=1e-9)
+    np.testing.assert_allclose(to_numpy(u.dagger()), expected.conj().T, atol=1e-9)
+    np.testing.assert_allclose(to_numpy(u + u), 2 * expected, atol=1e-9)
+    np.testing.assert_allclose(to_numpy(u + u - u), expected, atol=1e-9)
+    with pytest.raises(ValueError, match="different k"):
+        u + Matrix.x(FiveTuple)
+
+
 def test_sum_and_difference():
     np.testing.assert_allclose(to_numpy(Matrix.x() + Matrix.z()), X + Z)
     np.testing.assert_allclose(to_numpy(Matrix.x() - Matrix.z()), X - Z)
@@ -137,6 +166,15 @@ def test_density_from_vector():
     np.testing.assert_allclose(to_numpy(Matrix.density_from_vector(plus)), np.full((2, 2), 0.5), atol=1e-9)
 
 
+def test_density_from_vector_keeps_fivetuple_scale():
+    # #82, F59: |+> = (1, 1) / sqrt(2)^1, so |+><+| = (1, 1; 1, 1) / sqrt(2)^2
+    plus = Vector(q=2, element_representation=FiveTuple, k=1)
+    plus[0], plus[1] = FiveTuple.one(None), FiveTuple.one(None)
+    density = Matrix.density_from_vector(plus)
+    assert density.k == 2
+    np.testing.assert_allclose(to_numpy(density), np.full((2, 2), 0.5), atol=1e-9)
+
+
 def test_five_tuple_helpers():
     matrix = Matrix.x(FiveTuple).multiply_by_omega()
     np.testing.assert_allclose(to_numpy(matrix), OMEGA * X, atol=1e-9)
@@ -148,7 +186,7 @@ def test_five_tuple_helpers():
 
 
 def test_repr():
-    assert repr(Matrix.x()) == "[(0 + 0j), (1 + 0j), ]\n[(1 + 0j), (0 + 0j), ]\n"
+    assert repr(Matrix.x()) == "[(0 + 0j), (1 + 0j), ]\n[(1 + 0j), (0 + 0j), ]\nk=0"
 
 
 def test_matrix_of_variables():