From a98519b08eaf456dd74869d01233b3767575a142 Mon Sep 17 00:00:00 2001 From: Pete Shadbolt Date: Mon, 15 Feb 2016 15:50:13 +0000 Subject: [PATCH] Add edgelist method --- abp.py | 30 ++- graphsim.cpp | 528 ++++++++++++++++++++++++++++++++++++++++++++ tests/test_graph.py | 16 ++ 3 files changed, 567 insertions(+), 7 deletions(-) create mode 100644 graphsim.cpp create mode 100644 tests/test_graph.py diff --git a/abp.py b/abp.py index fa51621..3823a63 100644 --- a/abp.py +++ b/abp.py @@ -1,28 +1,44 @@ class Vertex(object): - def __init__(self): + + def __init__(self, index): + self.index = index self.neighbors = set() + def edgelist(self): + return [set((self.index, n)) for n in self.neighbors] + def __str__(self): - return "v: {}".format(", ".join(map(str, self.neighbors))) + return "{}".format(", ".join(map(str, self.neighbors))) + class GraphRegister(object): + def __init__(self, n): - self.vertices = [Vertex() for i in xrange(n)] + self.vertices = [Vertex(i) for i in xrange(n)] def add_edge(self, v1, v2): self.vertices[v1].neighbors.add(v2) self.vertices[v2].neighbors.add(v1) - def del_edge(self, v1, v2): self.vertices[v1].neighbors.remove(v2) self.vertices[v2].neighbors.remove(v1) + def toggle_edge(self, v1, v2): + if v2 in self.vertices[v1].neighbors: + self.del_edge(v1, v2) + else: + self.add_edge(v1, v2) + + def edgelist(self): + return map(tuple, frozenset(frozenset((v.index, n)) + for v in self.vertices + for n in v.neighbors)) + def __str__(self, ): - return "\n".join(str(v) for v in self.vertices) + return "\n".join(str(v) for v in self.vertices if len(v.neighbors) > 0) if __name__ == '__main__': g = GraphRegister(10) - g.add_edge(0,1) + g.toggle_edge(0, 1) print g - diff --git a/graphsim.cpp b/graphsim.cpp new file mode 100644 index 0000000..e99beff --- /dev/null +++ b/graphsim.cpp @@ -0,0 +1,528 @@ +// graphsim.cpp + +#include "graphsim.h" + +#include +#include + +//!Random coin toss. Change here to insert your favorite RNG. +int bool_rand (void) { + return random () > RAND_MAX/2; +} + + +//! Instantiate a quantum register with 'numQubits' qubits, initally all in state |0>. +/*! If randomize > -1 the RNG will be seeded with the current time plus the value +of randomize. (Otherwise, it is not seeded.) That the value of randomize is +added to the seed is useful in parallel processing settings where you want to +ensure different seeds. (If you call this from Python, remember, that Python's +RNG is not seeded.) */ +GraphRegister::GraphRegister (VertexIndex numQubits, int randomize) + : vertices (numQubits) +{ + if (randomize > -1) { + srandom (time(NULL) + randomize); + } +} + +//! Copy constructor +/*! Clones a register */ +GraphRegister::GraphRegister (GraphRegister& gr) + : vertices (gr.vertices) +{ +} + +//! Add an edge to the graph underlying the state. +void GraphRegister::add_edge (VertexIndex v1, VertexIndex v2) { + assert (v1 != v2); + vertices[v1].neighbors.insert (v2); + vertices[v2].neighbors.insert (v1); + //D cerr << "adding edge " << v1 << " - " << v2 << endl; +} + +//! Delete an edge to the graph underlying the state. +void GraphRegister::del_edge (VertexIndex v1, VertexIndex v2) { + DBGOUT ("deling edge " << v1 << " - " << v2 << endl); + vertices[v1].neighbors.erase (v2); + vertices[v2].neighbors.erase (v1); +} + +//! Toggle an edge to the graph underlying the state. +/*! (i.e. add it if not present, and delete it if present.) */ +void GraphRegister::toggle_edge (VertexIndex v1, VertexIndex v2) +{ + int n1 = vertices[v1].neighbors.erase (v2); + if (n1 == 1) { + vertices[v2].neighbors.erase (v1); + } else { + assert (n1 == 0); + add_edge (v1, v2); + } +} + +//! Create the Stabilizer of the state. +/*! This is useful to print out the stabilizer (or to compare with CHP). +You can also use print_stabilizer.*/ +Stabilizer& GraphRegister::get_full_stabilizer (void) const +{ + hash_set all_qubits; + for (VertexIterConst i = vertices.begin(); i != vertices.end(); i++) { + all_qubits.insert (i-vertices.begin()); + } + Stabilizer *s = new Stabilizer (*this, all_qubits); + return *s; +} + +//! Prints out the description of the current state +/*! in terms of adjacency lists of the graph and the VOps.*/ +void GraphRegister::print_adj_list (ostream& os) const +{ + for (VertexIndex i = 0; i < vertices.size(); i++) { + print_adj_list_line (os, i); + } +} + +//! Prints the line for Vertex i in the adjacency list representation of the state. +void GraphRegister::print_adj_list_line (ostream& os, VertexIndex i) const +{ + os << "Vertex " << i << ": VOp " + << vertices[i].byprod.get_name() << ", neighbors "; + for (VtxIdxIterConst j = vertices[i].neighbors.begin(); + j != vertices[i].neighbors.end(); j++) { + os << *j << " "; + } + os << endl; +} + +//! Print the current state in stabilizer representation. +void GraphRegister::print_stabilizer (ostream &os) const +{ + get_full_stabilizer().print (os); +} + +/*! Use the cphase look-up table. This is called by cphase after VOps that do not +commute with the cphase gate have been removed as far as possible. */ +void GraphRegister::cphase_with_table (VertexIndex v1, VertexIndex v2) +{ + // Structure of the table: + // first index: whether there was an edge between the operands before + // (0=no, 1= yes) + // second and third index: byprod op of v1 and v2 + // third index: information to obtain: + // 0= whether after the cphase there is an edges + // 1,2= new values of the byprod ops of v1 and v2Id + static const unsigned short cphase_tbl[2][24][24][3] = + #include "cphase.tbl" + ; + #ifdef DEBUGOUTPUT + DBGOUT ("cphase_with_table called on:\n"); + print_adj_list_line (cout, v1); + print_adj_list_line (cout, v2); + #endif + ConnectionInfo ci = getConnectionInfo (v1, v2); + unsigned op1 = vertices[v1].byprod.op; + unsigned op2 = vertices[v2].byprod.op; + // The table must only be used if a vertex has either no + // non-operand neighbors, or a diagonal byprod op + assert ((!ci.non1) || vertices[v1].byprod.is_diagonal()); + assert ((!ci.non2) || vertices[v2].byprod.is_diagonal()); + if (cphase_tbl[ci.wasEdge][op1][op2][0]) { + add_edge (v1, v2); + DBGOUT ("adding edge" << endl); + } else { + del_edge (v1, v2); + DBGOUT ("deling edge" << endl); + } + vertices[v1].byprod.op = cphase_tbl[ci.wasEdge][op1][op2][1]; + vertices[v2].byprod.op = cphase_tbl[ci.wasEdge][op1][op2][2]; + #ifdef DEBUGOUTPUT + DBGOUT ("cphase_with_table: after:\n"); + print_adj_list_line (cout, v1); + print_adj_list_line (cout, v2); + #endif + // The condition above must also hold afterwards: + ci = getConnectionInfo (v1, v2); + assert ((!ci.non1) || vertices[v1].byprod.is_diagonal()); + assert ((!ci.non2) || vertices[v2].byprod.is_diagonal()); +} + +/*! Check whether the qubits are connected to each other and to non-operand vertices.*/ +ConnectionInfo GraphRegister::getConnectionInfo + (VertexIndex v1, VertexIndex v2) { + ConnectionInfo ci; + ci.wasEdge = + vertices[v1].neighbors.find(v2) != vertices[v1].neighbors.end(); + if (! ci.wasEdge) { + ci.non1 = vertices[v1].neighbors.size() >= 1; + ci.non2 = vertices[v2].neighbors.size() >= 1; + } else { + ci.non1 = vertices[v1].neighbors.size() >= 2; + ci.non2 = vertices[v2].neighbors.size() >= 2; + } + return ci; +} + +//! Do a conditional phase gate between the two qubits. +void GraphRegister::cphase (VertexIndex v1, VertexIndex v2) +{ + // If there are non-operand neighbors, we can use neighborhood inversion + // to remove the byprod operators. + // These will store whether the operand vertices have nonoperand neighbors. + #ifdef DEBUGOUTPUT + DBGOUT ("before cphase between " << v1 << " and " << v2 << endl); + print_adj_list_line (cout, v1); + print_adj_list_line (cout, v2); + #endif + ConnectionInfo ci = getConnectionInfo (v1, v2); + + if (ci.non1) { + DBGOUT ("cphase: left vertex has NONs -> putting it to Id\n"); + remove_byprod_op (v1, v2); + } + ci = getConnectionInfo (v1, v2); + if (ci.non2) { + DBGOUT ("cphase: right vertex has NONs -> putting it to Id\n"); + remove_byprod_op (v2, v1); + } + ci = getConnectionInfo (v1, v2); + if (ci.non1 && !vertices[v1].byprod.is_diagonal()) { + // this can happen if v1 was first skipped + DBGOUT ("cphase: left one needs treatment again -> putting it to Id\n"); + remove_byprod_op (v1, v2); + } + cphase_with_table (v1, v2); +} + + +//! Do a controlled not gate between the vertices vc (control) and vt (target). +void GraphRegister::cnot (VertexIndex vc, VertexIndex vt) { + hadamard (vt); + cphase (vc, vt); + hadamard (vt); +} + +//! This structure is needed only by toggle_edges +struct Edge: public pair { + Edge (const VertexIndex a, const VertexIndex b) { + if (a < b) { + first = a; second = b; + } else { + first = b; second = a; + } + } +}; + +//! This structure is needed only by toggle_edges +struct edge_hash { + size_t operator() (const Edge& e) const + { + return e.first << 16 ^ e.second; + }; +}; + + +//! Toggles the edges between the vertex sets vs1 and vs2. +/*! The function takes extra care not to invert an edge twice. If vs1 and +vs2 are disjunct, this cannot happen and we do not need the function. +If vs1 == v2s, we can do without, too. */ +void GraphRegister::toggle_edges (const hash_set vs1, + const hash_set vs2) +{ + hash_set procd_edges; + for (VtxIdxIterConst i = vs1.begin(); i != vs1.end(); i++) { + for (VtxIdxIterConst j = vs2.begin(); j != vs2.end(); j++) { + if ((*i != *j) && + (procd_edges.find (Edge (*i, *j)) == procd_edges.end())) { + procd_edges.insert (Edge (*i, *j)); + toggle_edge (*i, *j); + } + } + } +} + +//! Get the adjacency matrix as an edgelist +std::vector GraphRegister::adj_list () +{ + std::vector out; + for (int i = 0; i < 100; ++i) + { + out.push_back(i); + } + return out; +} + +//! Measure the bare graph state in the Z basis. +int GraphRegister::graph_Z_measure (VertexIndex v, int force) +{ + int res; + if (force == -1) { + res = bool_rand (); + } else { + res = force; + } + DBGOUT ("gZm" << v << "," << res << " "); + #ifdef DEBUGOUTPUT + print_adj_list_line (cout, v); + #endif + hash_set nbg = vertices[v].neighbors; + for (VtxIdxIter i = nbg.begin(); i != nbg.end(); i++) { + del_edge (v, *i); + if (res) { + vertices[*i].byprod = vertices[*i].byprod * lco_Z; + } + } + if (! res) { + vertices[v].byprod = vertices[v].byprod * lco_H; + } else { + vertices[v].byprod = vertices[v].byprod * lco_X * lco_H; + } + return res; +} + +//! Measure the bare graph state in the Y basis. +int GraphRegister::graph_Y_measure (VertexIndex v, int force) +{ + int res; + if (force == -1) { + res = bool_rand (); + } else { + res = force; + } + DBGOUT ("gYm" << v << "," << res << " "); + hash_set vnbg = vertices[v].neighbors; + for (VtxIdxIter i = vnbg.begin(); i != vnbg.end(); i++) { + if (res) { + vertices[*i].byprod = vertices[*i].byprod * lco_spiZ; + } else { + vertices[*i].byprod = vertices[*i].byprod * lco_smiZ; + } + } + vnbg.insert (v); // Now, vnbg is the set of v and its neighbours. + for (VtxIdxIter i = vnbg.begin(); i != vnbg.end(); i++) { + for (VtxIdxIter j = i; j != vnbg.end(); j++) { + if (i != j) { + toggle_edge (*i, *j); + } + } + } + if (! res) { + // Messergebnis: +|0y> + vertices[v].byprod = vertices[v].byprod * lco_S; + } else { + // Messergebnis: -|0y> + vertices[v].byprod = vertices[v].byprod * lco_S.herm_adjoint(); + } + return res; +} + +//! Measure the bare graph state in the X basis. +int GraphRegister::graph_X_measure (VertexIndex v, bool* determined, + int force) +{ + if (vertices[v].neighbors.size () == 0) { + //not entangled qubit => result always 0 + DBGOUT ("gXm" << v << ",D "); + if (determined != NULL) + *determined = true; + return 0; + } + if (determined != NULL) + *determined = false; + // entangled qubit => let's get on with the complicated procedure + // throw a die: + int res; + if (force == -1) { + res = bool_rand (); + } else { + res = force; + } + DBGOUT ("gXm" << v << "," << res << " "); + DBGOUT ("forced: " << force << endl); + VertexIndex vb = *vertices[v].neighbors.begin(); // the choosen vertex + //D cerr << "vb = " << vb << endl; + // preparation step: store the neighborhood of v and vb + hash_set vn = vertices[v].neighbors; + hash_set vbn = vertices[vb].neighbors; + // First, put the byproduct ops: + if (! res) { + // measured a |+>: + // lco_spiY on vb + vertices[vb].byprod = vertices[vb].byprod * lco_spiY; + // Z on all in nbg(v) \ nbg(vb) \ {vb} + for (VtxIdxIter i = vertices[v].neighbors.begin(); + i != vertices[v].neighbors.end(); i++) { + if (*i != vb && + vertices[vb].neighbors.find(*i) == vertices[vb].neighbors.end()) { + vertices[*i].byprod = vertices[*i].byprod * lco_Z; + } + } + } else { + // measured a |->: + // lco_smiY on vb, and lco_Z on v: + vertices[vb].byprod = vertices[vb].byprod * lco_smiY; + vertices[v].byprod = vertices[v].byprod * lco_Z; + // Z on all in nbg(vb) \ nbg(v) \ {v} + for (VtxIdxIter i = vertices[vb].neighbors.begin(); + i != vertices[vb].neighbors.end(); i++) { + if (*i != v && + vertices[v].neighbors.find(*i) == vertices[v].neighbors.end()) { + vertices[*i].byprod = vertices[*i].byprod * lco_Z; + } + } + } + // Toggling the edges in three steps + // STEP 1: complement with Edges (nbg(v), nbg(vb)): + toggle_edges (vn, vbn); + // STEP 2: complement with the complete subgraph induced by the + // intersection of nbg(v) and nbg(vb): + // First, make the intersection + hash_set isc; + for (VtxIdxIter i = vn.begin(); i != vn.end(); i++) { + if (vbn.find(*i) != vbn.end()) { + isc.insert (*i); + } + } + // Now, toggle the edges + for (VtxIdxIter i = isc.begin(); i != isc.end(); i++) { + for (VtxIdxIter j = i; j != isc.end(); j++) { + if (*i != *j) { + toggle_edge (*i, *j); + } + } + } + // STEP 3: Toggle all edges from vb to nbg(v) \ {vb} + for (VtxIdxIter i = vn.begin(); i != vn.end(); i++) { + if (*i != vb) { + toggle_edge (vb, *i); + } + } + return res; +} + +//! Measure qubit v in basis 'basis'. +/* For basis, pass a LocCliffOp object which has to be equal to one of lco_X, lco_Y +or lco_Z. If you want to now, whether the result was choosen at random or determined +by the state, pass a bool pointer in which this information will be written. +If you want to force the result to be a certain value, pass 0 or 1 to 'force'. This only +works, if the result is not determined. If it is, 'force' is ignored.*/ +int GraphRegister::measure (VertexIndex v, LocCliffOp basis, bool* determined, + int force) +{ + assert (basis.op >= lco_X.op && basis.op <= lco_Z.op); + if (determined != NULL) { + *determined = false; + } + assert (force >= -1 && force <= 1); + LocCliffOp basis_orig = basis; + RightPhase rp = basis.conjugate (vertices[v].byprod.herm_adjoint()); + assert (rp == rp_p1 || rp == rp_m1); + if (force != -1 && rp == rp_m1) { + force = force ^ 0x01; + } + int res; + switch (basis.op) { + case 1 /* lco_X */: res = graph_X_measure (v, determined, force); break; + case 2 /* lco_Y */: res = graph_Y_measure (v, force); break; + case 3 /* lco_Z */: res = graph_Z_measure (v, force); break; + default: exit (1); + } + if (rp == rp_m1) { + res = ! res; + } else { + assert (rp == rp_p1); + } + // check: the measured vertex should be singled out: + assert (vertices[v].neighbors.size() == 0); + // Check that the vertex is now in the correct eigenstate: + LocCliffOp assert_op = lco_X; + assert (assert_op.conjugate (vertices[v].byprod) == (res ? rp_m1 : rp_p1)); + assert (assert_op == basis_orig); + return res; +} + +//! Do a neighborhood inversion (i.e. local complementation) about vertex v. +/* This changes the state's graph representation but not the state itself, as the +necessary correction to the VOps are applied. */ +void GraphRegister::invert_neighborhood (VertexIndex v) +{ + // Invert the neighborhood: + #ifdef DEBUGOUTPUT + DBGOUT ("Inverting about "); + print_adj_list_line (cout, v); + #endif + hash_set vn = vertices[v].neighbors; + for (VtxIdxIter i = vn.begin(); i != vn.end(); i++) { + for (VtxIdxIter j = i; j != vn.end(); j++) { + if (*i != *j) { + //cerr << "toggling " << *i << "," << *j << endl; + toggle_edge (*i, *j); + } + } + // and adjust the local Cliffords: + vertices[*i].byprod = vertices[*i].byprod * lco_spiZ.herm_adjoint(); + } + // finally, adjust the local Clifford of v: + vertices[v].byprod = vertices[v].byprod * lco_smiX.herm_adjoint(); +} + +//! Do neighborhood inversions to reduce the VOp of vertex v to the identity. +/* 'avoid' is avoided as swapping partner, i.e. the swapping partner will not +be 'avoid' unless this is the only neighbor of v. If no neighbor is available, +the program exits with error message.*/ +bool GraphRegister::remove_byprod_op (VertexIndex v, VertexIndex avoid) +{ + // A lookup table on how any LC operator can be composed from them + // generators lco_spiZ (denoted 'V') and lco_smiX (denoted 'U'): + // (The lookup table was generated from sqrt-gen.nb and copied here manually.) + const char * comp_tbl [24] = + {"UUUU", "UU", "VVUU", "VV", + "VUU", "V", "VVV", "UUV", + "UVU", "UVUUU", "UVVVU", "UUUVU", + "UVV", "VVU", "UUU", "U", + "VVVU", "UUVU", "VU", "VUUU", + "UUUV", "UVVV", "UV", "UVUU"}; + // Of course, we need a neighborhood + if (vertices[v].neighbors.size() == 0) { + cerr << "remove_byprod_op: called with isolated vertex. Aborting.\n"; + } + // This will be the swapping partner: + VertexIndex vb = * vertices[v].neighbors.begin (); + if (vb == avoid) { + // Is there an alternative to 'avoid'? If so, use it. + VtxIdxIter vbi = vertices[v].neighbors.begin (); + vbi++; + if (vbi != vertices[v].neighbors.end()) { + vb = *vbi; + } + } + #ifdef DEBUGOUTPUT + DBGOUT ("remove_byprod_op called: (v, avoid, vb):\n"); + print_adj_list_line (cout, v); + print_adj_list_line (cout, avoid); + print_adj_list_line (cout, vb); + #endif + const char * comp = comp_tbl[vertices[v].byprod.op]; + DBGOUT ("using " << comp << endl); + for (int pos = strlen(comp)-1; pos >= 0; pos--) { + if (comp[pos] == 'U') { + // A U will vanish if we do an inversion on v + DBGOUT ("U ->"); + invert_neighborhood (v); + } else { + assert (comp[pos] == 'V'); + DBGOUT ("V ->"); + // For this we need to invert on a neighbor of v + invert_neighborhood (vb); + } + } + #ifdef DEBUGOUTPUT + DBGOUT ("remove_byprod_op, after: (v, avoid, vb):\n"); + print_adj_list_line (cout, v); + print_adj_list_line (cout, avoid); + print_adj_list_line (cout, vb); + #endif + // Now, we should have lco_Id left + assert (vertices[v].byprod == lco_Id); + return true; +} + diff --git a/tests/test_graph.py b/tests/test_graph.py new file mode 100644 index 0000000..cb4ef24 --- /dev/null +++ b/tests/test_graph.py @@ -0,0 +1,16 @@ +from nose import with_setup +import abp + +def setup(): + global g + g = abp.GraphRegister(10) + +@with_setup(setup) +def test_adding(): + g.add_edge(0,1) + g.add_edge(1,2) + g.toggle_edge(2,0) + g.toggle_edge(2,3) + g.toggle_edge(3,2) + print list(g.edgelist()) +