Merge branch 'measure'

complex.c

@@ -35,6 +35,10 @@ int quda_complex_eq(complex_t op1, complex_t op2) {
 	return 0;
 }
 
+float quda_complex_abs_square(complex_t c) {
+	return c.real*c.real + c.imag*c.imag;
+}
+
 float quda_complex_abs(complex_t c) {
 	float res = c.real*c.real + c.imag*c.imag;
 	return sqrt(res);

complex.h

@@ -32,6 +32,13 @@ complex_t quda_complex_copy(complex_t c);
 /* Test equality of two complex numbers (returns 1 if equal, 0 otherwise) */
 int quda_complex_eq(complex_t op1, complex_t op2);
 
+/* Complex absolute square
+ * Equivalent to a complex number times its conjugate
+ * For a quantum state amplitude, this value represents the probability
+ * of the quantum state.
+ */
+float quda_complex_abs_square(complex_t c);
+
 /* Complex modulus (absolute value) */
 float quda_complex_abs(complex_t c);
 

quantum_reg.c

@@ -3,6 +3,7 @@
 
 #include "quantum_reg.h"
 #include <stdlib.h>
+#include <math.h>
 
 int quda_quantum_reg_init(quantum_reg* qreg, int qubits) {
 	qreg->qubits = qubits;
@@ -26,6 +27,125 @@ void quda_quantum_reg_delete(quantum_reg* qreg) {
 	free(qreg->states);
 }
 
+void quda_quantum_bit_set(int target, quantum_reg* qreg) {
+	int i;
+	uint64_t mask = 1 << target;
+	for(i=0;i<qreg->num_states;i++) {
+		qreg->states[i].state = qreg->states[i].state | mask;
+	}
+
+	quda_quantum_reg_coalesce(qreg);
+}
+
+void quda_quantum_bit_reset(int target, quantum_reg* qreg) {
+	int i;
+	uint64_t mask = ~(1 << target);
+	for(i=0;i<qreg->num_states;i++) {
+		qreg->states[i].state = qreg->states[i].state & mask;
+	}
+
+	quda_quantum_reg_coalesce(qreg);
+}
+
+/* Performs a measurement on the quantum register */
+int quda_quantum_reg_measure(quantum_reg* qreg, uint64_t* retval) {
+	if(retval == NULL) return -2;
+	float f = quda_rand_float();
+	int i;
+	for(i=0;i<qreg->num_states;i++) {
+		if(!quda_complex_eq(qreg->states[i].amplitude,QUDA_COMPLEX_ZERO)) {
+			f -= quda_complex_abs_square(qreg->states[i].amplitude);
+			if(f < 0) {
+				*retval = qreg->states[i].state;
+				return 0;
+			}
+		}
+	}
+
+	return -1;
+}
+
+/* Performs a real-world quantum measurement.
+ * The register collapses to the physical state measured with probability 1.
+ */
+int quda_quantum_reg_measure_and_collapse(quantum_reg* qreg, uint64_t* retval) {
+	if(retval == NULL) return -2;
+	float f = quda_rand_float();
+	int i;
+	for(i=0;i<qreg->num_states;i++) {
+		if(!quda_complex_eq(qreg->states[i].amplitude,QUDA_COMPLEX_ZERO)) {
+			f -= quda_complex_abs_square(qreg->states[i].amplitude);
+			if(f < 0) {
+				*retval = qreg->states[i].state;
+				qreg->states[0].state = qreg->states[i].state;
+				qreg->states[0].amplitude = QUDA_COMPLEX_ONE;
+				qreg->num_states = 1;
+				return 0;
+			}
+		}
+	}
+
+	return -1;
+}
+
+/* Measure 1 bit of a quantum register */
+int quda_quantum_bit_measure(int target, quantum_reg* qreg) {
+	float p = 0;
+	float f = quda_rand_float();
+	uint64_t mask = 1 << target;
+	int i;
+	// Accumulate probability that the bit is in state |1>
+	for(i = 0;i<qreg->num_states;i++) {
+		if(qreg->states[i].state & mask) {
+			p += quda_complex_abs_square(qreg->states[i].amplitude);
+			// TODO: Determine overhead of this comparison
+			if(p > quda_rand_float()) return 1;
+		}
+	}
+
+	return 0;
+	//return quda_rand_float() < p ? 1 : 0;
+}
+
+int quda_quantum_bit_measure_and_collapse(int target, quantum_reg* qreg) {
+	// Measure bit conventionally
+	// TODO: Can allow probability measurement from conventional to complete and remove it below
+	int retval = quda_quantum_bit_measure(target,qreg);
+
+	// Collapse states to those possible
+	uint64_t mask = 1 << target;
+	float p = 0;
+	int i;
+	for(i=0;i<qreg->num_states;i++) {
+		// TODO: Ideally, remove nested conditions
+		// TODO: Actually prune in this loop instead of just invalidating
+		if(qreg->states[i].state & mask) {
+			if(retval) { // this is a valid state, accumulate probability to renormalize
+				p += quda_complex_abs_square(qreg->states[i].amplitude);
+			} else { // this is an invalid state -- nullify
+				qreg->states[i].amplitude = QUDA_COMPLEX_ZERO;
+			}
+		} else {
+			if(!retval) { // valid state, accumulate probability
+				p += quda_complex_abs_square(qreg->states[i].amplitude);
+			} else { // invalid state -- nullify
+				qreg->states[i].amplitude = QUDA_COMPLEX_ZERO;
+			}
+		}
+	}
+
+	// TODO: Remove this call for optimization within the above loop
+	quda_quantum_reg_prune(qreg);
+
+	// Renormalize
+	float k = sqrt(1.0f/p);
+	for(i=0;i<qreg->num_states;i++) {
+		qreg->states[i].amplitude = quda_complex_rmul(qreg->states[i].amplitude,k);
+	}
+
+	return retval;
+}
+
 void quda_quantum_reg_prune(quantum_reg* qreg) {
 	int i,end;
 	for(i=0,end=qreg->num_states-1;i < end;i++) {
@@ -112,6 +232,10 @@ int quda_quantum_reg_trim(quantum_reg* qreg) {
 	return 0;
 }
 
+float quda_rand_float() {
+	return rand()/(float)RAND_MAX;
+}
+
 int qstate_compare(const void* qstate1, const void* qstate2) {
 	uint64_t diff = ((quantum_state_t*)qstate1)->state - ((quantum_state_t*)qstate2)->state;
 	if(diff == 0) return 0;

quantum_reg.h

@@ -29,14 +29,54 @@ typedef struct quantum_reg {
  */
 int quda_quantum_reg_init(quantum_reg* qreg, int qubits);
 
-/* Sets the register to a single physical state with probability 1. */
-void quda_quantum_reg_set(quantum_reg* qreg, uint64_t state);
-
 /* Frees the given register for deletion. 
  * If the passed qreg was dynamically allocated, it must still be freed separately.
  */
 void quda_quantum_reg_delete(quantum_reg* qreg);
 
+/* Sets the register to a single physical state with probability 1. */
+void quda_quantum_reg_set(quantum_reg* qreg, uint64_t state);
+
+/* Sets a single bit of a quantum register to 1 with probability 1.
+ * If the system is in a superposition, it is collapsed into the subset
+ * of possible states allowed by this value.
+ */
+void quda_quantum_bit_set(int target, quantum_reg* qreg);
+
+/* Sets a single bit of a quantum register to 0 with probability 1.
+ * If the system is in a superposition, it is collapsed into the subset
+ * of possible states allowed by this value.
+ */
+void quda_quantum_bit_reset(int target, quantum_reg* qreg);
+
+/* Performs a measurement on the quantum register and stores the state
+ * in 'retval' if non-NULL.
+ * Returns 0 on success, -1 on retval NULL, -2 on normalization error.
+ */
+// TODO: Attempt correction for minor normalization errors (ie floating point precision errors)
+int quda_quantum_reg_measure(quantum_reg* qreg, uint64_t* retval);
+
+/* Performs a real-world quantum measurement and stores the state in 
+ * 'retval' if non-NULL.
+ * On success, returns 0 and the register collapses to the physical state
+ * measured with probability 1.
+ * On failure, does not collapse. Returns -1 on retval NULL, -2 on
+ * normalization error.
+ */
+// TODO: Attempt correction for minor normalization errors (ie floating point precision errors)
+int quda_quantum_reg_measure_and_collapse(quantum_reg* qreg, uint64_t* retval);
+
+/* Measure 1 bit of a quantum register */
+int quda_quantum_bit_measure(int target, quantum_reg* qreg);
+
+/* Perform a real-world quantum measurement of 1 quantum register bit.
+ * The system collapses into the subset of possible states allowed by
+ * the value of the measured bit.
+ * Simultaneously prunes zero-amplitude states.
+ * Does not coalesce identical states.
+ */
+int quda_quantum_bit_measure_and_collapse(int target, quantum_reg* qreg);
+
 /* Removes zero-amplitude states from the register. */
 void quda_quantum_reg_prune(quantum_reg* qreg);
 
@@ -59,6 +99,11 @@ void quda_quantum_reg_coalesce(quantum_reg* qreg);
  */
 int quda_quantum_reg_trim(quantum_reg* qreg);
 
+/* Generates a float in the range [0,1) */
+// TODO: Look at performance implications of using 'double' here
+float quda_rand_float();
+
+/* Comparator for sorting the quantum_state array */
 int qstate_compare(const void* qstate1, const void* qstate2);
 
 #endif // __QUDA_QUANTUM_REG_H