Updated validation.py

qaoa/utils/validation.py

@@ -9,17 +9,27 @@ def _bitstring(i, n, flip=False):
     else:
         return format(i, f'0{n}b')
 
+def phase_cost(problem, bitstring, global_phase, omit_invalid_states, flip):
+    idx = int(bitstring[::-1] if flip else bitstring, 2)
+    if omit_invalid_states and not problem.isFeasible(bitstring):
+        energies[idx] = global_phase
+        mask[idx] = False
+        return
+    energies[idx] = problem.energy(bitstring)
 
 def check_phase_separator_exact_qaoa(qaoa, *arg, **kwarg):
     return check_phase_separator_exact_problem(qaoa.problem, *arg, **kwarg)
 
-def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-8, rtol=1e-8):
+def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-7, rtol=1e-7, global_phase = 0.0, omit_invalid_states = False):
     """
-    Exact check that the problem's circuit represents the problem's energy function.
+    Exact check that the problem's circuit represents the problem's cost function.
     This tests checks that the unitary operator represented by the quantum circuit is
     equal to the expected matrix with diagonal elements 
-    exp(-j*t*energy(e)),
+    exp(-j*t*cost(e)),
     where e is the corresponding binary state, up to a global phase.
+
+    Works for anzats' that remain in a feasible subspace. If the anzats' can reach infeasible spaces, 
+    the problem must define a suitable penalty for such solutions.
 
     Suitable for <= 10 qubits as this check uses the full unitary matrix of size 2^n x 2^n).
     Returns: (ok: bool, report: dict)
@@ -30,31 +40,33 @@ def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-8, rtol
         {problem.circuit.parameters[0]: t},
         inplace = False
     )
-    energy_fn = problem.energy
 
     U = Operator(circ).data  # complex ndarray
+    diag = np.diag(U)
     n = circ.num_qubits
     d = 2**n
+
     # Compare diagonal phases to expected, modulo a global phase
     # expected diag entries
-    energies = []
+    global energies, mask
+
+    energies = [0]*len(diag)
+    mask = np.ones(len(diag), dtype=bool)
+
     for i in range(d):
-        energies.append(energy_fn(_bitstring(i, n, flip=flip)))
+        phase_cost(problem, _bitstring(i, n, flip=flip), global_phase, omit_invalid_states, flip=flip)
+
     expected = np.exp(-1j * t * np.asarray(energies, dtype=float))
-
 
-    diag = np.diag(U)
-    # if n < 4: 
-    #     for i in range(d):
-    #         print(expected[i]* diag[0], diag[i])
     # Remove global phase by aligning first nonzero expected
     ref_idx = 0
     g = diag[ref_idx] / expected[ref_idx]  # global phase factor
     ratios = diag / (expected * g)
 
     # Errors
-    mag_err = np.max(np.abs(np.abs(diag) - 1.0))
-    phase_err = np.max(np.abs(np.angle(ratios)))  # max residual phase after removing global
+    mag_err = np.max(np.abs(np.abs(diag[mask]) - 1.0))
+    phase_err = np.max(np.abs(np.angle(ratios[mask])))
+
     ok = (mag_err <= rtol) and (phase_err <= atol)
 
     report = {
@@ -63,6 +75,8 @@ def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-8, rtol
         "max_phase_error_rad_after_global": float(phase_err),
         "global_phase_rad": float(np.angle(g)),
     }
+
+
     if not ok:
         # include a few worst offenders
         idx_sorted = np.argsort(-np.abs(np.angle(ratios)))
@@ -71,9 +85,10 @@ def check_phase_separator_exact_problem(problem, t=1, flip=True, atol=1e-8, rtol
             bad.append({
                 "bitstring": list(_bitstring(k, n, flip=flip)),
                 "diag_entry": complex(diag[k]),
+                "raw expected": complex(expected[k]),
                 "expected": complex(expected[k]*g),
                 "phase_residual_rad": float(np.angle(ratios[k])),
                 "magnitude": float(np.abs(diag[k])),
             })
         report["examples"] = bad
-    return ok, report
+    return ok, report