From 8388c110e5dc01f18e85966c0c15a53a31dcaeca Mon Sep 17 00:00:00 2001 From: Sebastiaan Brand Date: Mon, 22 Jul 2024 14:33:41 +0200 Subject: [PATCH] added modifiable threshold to circuit equivalence check --- examples/circuit_equivalence.c | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/examples/circuit_equivalence.c b/examples/circuit_equivalence.c index 7030d22..c8ba759 100644 --- a/examples/circuit_equivalence.c +++ b/examples/circuit_equivalence.c @@ -8,6 +8,7 @@ #include "../qasm/qsylvan_qasm_parser.h" #define max(a,b) ((a > b) ? a : b) +#define min(a,b) ((a < b) ? a : b) /*********************************************/ @@ -19,6 +20,7 @@ static size_t min_cachesize = 1LL<<16; static size_t max_cachesize = 1LL<<18; static size_t wgt_tab_size = 1LL<<23; static double tolerance = 1e-14; +static double threshold = 1e-6; static int wgt_table_type = COMP_HASHMAP; static int wgt_norm_strat = NORM_MAX; static bool count_nodes = false; @@ -30,6 +32,7 @@ static struct argp_option options[] = {"workers", 'w', "", 0, "Number of workers/threads (default=1)", 0}, {"norm-strat", 's', "", 0, "Edge weight normalization strategy (default=max)", 0}, {"tol", 't', "", 0, "Tolerance for deciding edge weights equal (default=1e-14)", 0}, + {"threshold", 1000, "", 0, "Threshold for deciding when two matrices are close enough to be equivalent (default=1e-6)", 0}, {"count-nodes", 'c', 0, 0, "Track maximum number of nodes", 0}, {0, 0, 0, 0, 0, 0} }; @@ -53,6 +56,9 @@ parse_opt(int key, char *arg, struct argp_state *state) case 'c': count_nodes = true; break; + case 1000: + threshold = atof(arg); + break; case ARGP_KEY_ARG: if (state->arg_num == 0) circuit_U = parse_qasm_file(arg); @@ -77,6 +83,7 @@ static struct argp argp = { options, parse_opt, " ", 0, 0, 0, 0 typedef struct stats_s { bool equivalent; char counterexample[100]; + double fidelity; double cpu_time; double wall_time; size_t max_nodes_total; @@ -94,6 +101,7 @@ void print_stats() { printf(" \"circuit_V\": \"%s\",\n", circuit_V->name); printf(" \"counterexample\": \"%s\",\n", stats.counterexample); printf(" \"equivalent\" : %d,\n", (int)stats.equivalent); + printf(" \"min_fidelity\" : %.20lf\n", stats.fidelity); printf(" \"max_nodes_total\": %" PRIu64 ",\n", stats.max_nodes_total); printf(" \"max_nodes_U\": %" PRIu64 ",\n", stats.max_nodes_U); printf(" \"max_nodes_V\": %" PRIu64 ",\n", stats.max_nodes_V); @@ -219,6 +227,7 @@ QMDD compute_UPUdag(quantum_circuit_t *circuit, gate_id_t P, BDDVAR k) { BDDVAR nqubits = U->qreg_size; stats.equivalent = true; + stats.fidelity = 1.0; uint64_t m = 0; uint32_t XZ[2] = {GATEID_X, GATEID_Z}; @@ -229,7 +238,6 @@ QMDD compute_UPUdag(quantum_circuit_t *circuit, gate_id_t P, BDDVAR k) { aadd_protect(&qmdd_U); QMDD qmdd_V = compute_UPUdag(V, XZ[i], k); aadd_unprotect(&qmdd_U); - // TODO: ^ remove global phase? if (count_nodes) { size_t nodes_U = aadd_countnodes(qmdd_U); @@ -240,10 +248,18 @@ QMDD compute_UPUdag(quantum_circuit_t *circuit, gate_id_t P, BDDVAR k) { } if (qmdd_U != qmdd_V) { - stats.equivalent = false; - snprintf(stats.counterexample, sizeof(stats.counterexample), - "U*%c_%d*U^dag != V*%c_%d*V^dag", XZ_str[i], k, XZ_str[i], k); - return; + // if not exactly equal, compute overlap + double norm = pow(2.0, nqubits*2); + double fid = qmdd_fidelity(qmdd_U, qmdd_V, nqubits*2) / norm; + stats.fidelity = min(stats.fidelity, fid); + //printf("F(U*%c_%d*U^dag,V*%c_%d*V^dag) = %.10lf\n", XZ_str[i], k, XZ_str[i], k, fid); + + if (fabs(fid - 1.0) > threshold) { + stats.equivalent = false; + snprintf(stats.counterexample, sizeof(stats.counterexample), + "U*%c_%d*U^dag != V*%c_%d*V^dag", XZ_str[i], k, XZ_str[i], k); + return; + } } } }