Add nice hhl example

algorithms/quantum_linear_solvers/hhl/hhl.qmod

@@ -1,164 +1,110 @@
-qfunc load_b_expanded___0(amplitudes: real[], output memory: qbit[2]) {
-  prepare_amplitudes(amplitudes, 0.0, memory);
-}
+// Harrow–Hassidim–Lloyd (HHL) Algorithm for Solving Linear Systems
 
-qfunc apply_to_all_expanded___0(target: qbit[4]) {
-  repeat (index: 4) {
-    H(target[index]);
-  }
-}
+// This example implements the HHL quantum linear-systems algorithm to solve
+//   A |x⟩ = |b⟩
+// where A is a Hermitian matrix and |b⟩ is the right-hand-side vector prepared
+// as a quantum state.
 
-qfunc hamiltonian_evolution_with_power_0_lambda___0_0_expanded___0(pw: int, target: qbit[2]) {
-  power (pw) {
-    unitary([
-      [
-        ((-0.09406240950199857) + 0.8149069223122054j),
-        (0.03521871946675126 - 0.029763534641642615j),
-        ((-0.018800717000078293) - 0.16142879795007106j),
-        (0.43769245930764733 + 0.32705554908759304j)
-      ],
-      [
-        (0.03521871946675127 - 0.029763534641642633j),
-        ((-0.15347248298890326) - 0.17275282472948233j),
-        (0.23117644455908531 + 0.8872069971297389j),
-        (0.23971825754883572 + 0.21548267921288933j)
-      ],
-      [
-        ((-0.01880071700007826) - 0.16142879795007103j),
-        (0.23117644455908523 + 0.8872069971297386j),
-        ((-0.1219131720516462) + 0.13200138126428373j),
-        (0.29584069101495575 + 0.11488938733473114j)
-      ],
-      [
-        (0.43769245930764744 + 0.32705554908759327j),
-        (0.2397182575488357 + 0.21548267921288933j),
-        (0.29584069101495586 + 0.1148893873347311j),
-        ((-0.6563827949579104) + 0.25690988991104674j)
-      ]
-    ], target);
-  }
-}
 
-qfunc unitary_with_power_0_lambda___0_0_expanded___0(k: int, memory_captured__hhl__1: qbit[2]) {
-  hamiltonian_evolution_with_power_0_lambda___0_0_expanded___0(k, memory_captured__hhl__1);
-}
+// In this example, the matrix A is given as a sum of Pauli strings:
+//   A = 0.03*X(0)*X(1) + 0.05*Z(0)*Z(1) + 0.02*Y(0)*Y(1) + 0.08*I
+// This Pauli decomposition corresponds to the matrix
+//   A = [[0.13, 0.00, 0.00, 0.01],
+//        [0.00, 0.03, 0.05, 0.00],
+//        [0.00, 0.05, 0.03, 0.00],
+//        [0.01, 0.00, 0.00, 0.13]]
 
-qfunc qft_no_swap_expanded___0(qbv: qbit[4]) {
-  repeat (i: 4) {
-    H(qbv[i]);
-    repeat (j: (4 - i) - 1) {
-      CPHASE(pi / (2 ** (j + 1)), qbv[(i + j) + 1], qbv[i]);
-    }
-  }
-}
+hamiltonian: SparsePauliOp = SparsePauliOp {
+    terms=[
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=1, index=0},
+          IndexedPauli {pauli=1, index=1}
+        ],
+        coefficient=0.03
+      },
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=3, index=1},
+          IndexedPauli {pauli=3, index=0}
+        ],
+        coefficient=0.05
+      },
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=2, index=0},
+          IndexedPauli {pauli=2, index=1}
+        ],
+        coefficient=0.02
+      },
+      SparsePauliTerm {
+        paulis=[
+          IndexedPauli {pauli=0, index=0}
+        ],
+        coefficient=0.08
+      }
+    ],
+    num_qubits=2
+    };
 
-qfunc qft_expanded___0(target: qbit[4]) {
-  repeat (index: 2.0) {
-    SWAP(target[index], target[3 - index]);
-  }
-  qft_no_swap_expanded___0(target);
-}
+// For the right-hand-side vector |b⟩, we use a normalized vector of size 4.
+rhs_vector: real[] = [sqrt(1/6), sqrt(1/6), sqrt(1/6), sqrt(1/2)];
 
-qfunc qpe_flexible_expanded___0(phase: qbit[4], memory_captured__hhl__1: qbit[2]) {
-  apply_to_all_expanded___0(phase);
-  repeat (index: 4) {
-    control (phase[index]) {
-      unitary_with_power_0_lambda___0_0_expanded___0(2 ** index, memory_captured__hhl__1);
-    }
-  }
-  invert {
-    qft_expanded___0(phase);
-  }
+// In our example the matrix A is a sum of commuting Pauli strings.
+// We define how to take powers of the Hamiltonian evolution exp(2π A), which is applied
+// within the Quantum Phase Estimation. We simply multiply the evolution time by the power value,
+// keeping the number of Trotter repetition being 1, as suitable for commuting terms.
+qfunc powered_suzuki_trotter_commuting_terms(k: int, hamiltonian: SparsePauliOp, qba: qbit[]) {
+  suzuki_trotter(hamiltonian, -2*pi* k, 1, 1, qba);
 }
 
-qfunc assign_amplitude_table_expanded___0(const index: qbit[4], indicator: qbit) {
-  RY(0.49069646327199606, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.159660988575539, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.215103068840025, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.114786680603947, indicator);
-  skip_control {
-    CX(index[2], indicator);
-  }
-  RY(0.0734201414091853, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.222322935969005, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.181317761599329, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.22482930086683403, indicator);
-  skip_control {
-    CX(index[3], indicator);
-  }
-  RY(0.192520246080629, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.184296564729869, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.22312206827512002, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.0676093870745949, indicator);
-  skip_control {
-    CX(index[2], indicator);
-  }
-  RY(0.0978641117835104, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(-0.216731023385388, indicator);
-  skip_control {
-    CX(index[1], indicator);
-  }
-  RY(-0.168241915420623, indicator);
-  skip_control {
-    CX(index[0], indicator);
-  }
-  RY(0.309069995704199, indicator);
-  skip_control {
-    CX(index[3], indicator);
-  }
+// We define a quantum function that assigns eigenvalue inversion to the amplitudes for a signed qnum of
+// size 4 with 4 fraction places index. The desired amplitudes are conditioned by an additional
+// indicator qubit being at state 1, and are normalized by the smallest possible value of 1/2^index.size.
+qfunc assign_inversion_size4(const index: qbit[4], indicator: qbit) {
+  RY(-0.0157, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(-0.0157, indicator);
+  skip_control { CX(index[1], indicator); }
+  RY(-0.3379, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(0.3066, indicator);
+  skip_control { CX(index[2], indicator); }
+  RY(-0.1047, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(-0.1047, indicator);
+  skip_control { CX(index[1], indicator); }
+  RY(-0.3181, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(0.4649, indicator);
+  skip_control { CX(index[3], indicator); }
+  RY(-0.0475, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(-0.0475, indicator);
+  skip_control { CX(index[1], indicator); }
+  RY(-0.3407, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(0.2457, indicator);
+  skip_control { CX(index[2], indicator); }
+  RY(-0.0939, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(-0.0939, indicator);
+  skip_control { CX(index[1], indicator); }
+  RY(-0.3122, indicator);
+  skip_control { CX(index[0], indicator); }
+  RY(0.8154, indicator);
+  skip_control { CX(index[3], indicator); }
 }
 
-qfunc hhl_expanded___0(rhs_vector: real[], output memory: qbit[2], output estimator: qnum<4, False, 4>, output indicator: qbit) {
-  allocate(4, False, 4, estimator);
-  load_b_expanded___0([
-    0.18257418583505536,
-    0.3651483716701107,
-    0.7302967433402214,
-    0.5477225575051661
-  ], memory);
-  allocate(1, indicator);
+qfunc main(output solution: qnum<2, False, 0>, output phase_var: qnum<4, True, 4>, output indicator: qbit) {
+  allocate(phase_var);
+  allocate(indicator);
+  prepare_amplitudes(rhs_vector, 0.0, solution);
   within {
-    qpe_flexible_expanded___0(estimator, memory);
+    qpe_flexible(lambda(k) {
+      powered_suzuki_trotter_commuting_terms(k, hamiltonian, solution);
+    }, phase_var);
   } apply {
-    assign_amplitude_table_expanded___0(estimator, indicator);
+    assign_inversion_size4(phase_var, indicator);
   }
 }
-
-qfunc main(output res: qnum<2, False, 0>, output estimator_var: qnum<4, False, 4>, output indicator: qbit) {
-  hhl_expanded___0([
-    0.18257418583505536,
-    0.3651483716701107,
-    0.7302967433402214,
-    0.5477225575051661
-  ], res, estimator_var, indicator);
-}