Improved measurements implementation

aromanro · Feb 22, 2025 · 65951ac · 65951ac
1 parent 32a19ec
commit 65951ac
Show file tree

Hide file tree

Showing 4 changed files with 73 additions and 101 deletions.
diff --git a/QCSim/CheckCHSHInequality.h b/QCSim/CheckCHSHInequality.h
@@ -59,7 +59,7 @@ namespace BellInequalities {
 			int res1 = -1;
 			if (separateMeasurements)
 			{
-				state = BaseClass::Measure(0, 0);
+				state = BaseClass::Measure(0);
 				if (state) res1 = 1;
 			}
 
@@ -70,7 +70,7 @@ namespace BellInequalities {
 
 			int res2 = -1;
 			if (separateMeasurements)
-				state |= BaseClass::Measure(1, 1) << 1;
+				state |= BaseClass::Measure(1) << 1;
 			else
 			{
 				state = BaseClass::Measure();

diff --git a/QCSim/MeasurementsTests.cpp b/QCSim/MeasurementsTests.cpp
@@ -66,7 +66,7 @@ bool checkQubit0()
 
 		// measure first qubit only
 
-		const size_t state = reg.Measure(0, 0);
+		const size_t state = reg.MeasureQubit(0);
 		++measurements[state];
 
 		Eigen::VectorXcd amplitudes = reg.getRegisterStorage();
@@ -129,7 +129,7 @@ bool checkQubit1()
 
 		// measure second qubit only
 
-		const size_t state = reg.Measure(1, 1);
+		const size_t state = reg.MeasureQubit(1);
 		++measurements[state];
 
 		Eigen::VectorXcd amplitudes = reg.getRegisterStorage();
@@ -190,10 +190,10 @@ bool checkSingleQubitMeasurements()
 	{
 		setRegister(reg);
 
-		size_t state = reg.Measure(0, 0);
+		size_t state = reg.MeasureQubit(0);
 		++measurements[state];
 
-		state = reg.Measure(1, 1);
+		state = reg.MeasureQubit(1);
 		++fmeasurements[state];
 	}
 
@@ -211,10 +211,10 @@ bool checkSingleQubitMeasurements()
 	{
 		setRegister(reg);
 
-		size_t state = reg.Measure(1, 1);
+		size_t state = reg.MeasureQubit(1);
 		++measurements[state];
 
-		state = reg.Measure(0, 0);
+		state = reg.MeasureQubit(0);
 		++fmeasurements[state];
 	}
 

diff --git a/QCSim/QubitRegisterCalculator.h b/QCSim/QubitRegisterCalculator.h
@@ -299,15 +299,20 @@ namespace QC {
 
 			const size_t measuredQubitMask = 1ULL << qubit;
 
-			for (size_t state = measuredQubitMask; state < NrBasisStates; ++state)
+			size_t state = 0;
+			for (size_t i = 0; i < NrBasisStates; ++i)
 			{
-				if (state & measuredQubitMask)
-					accum += std::norm(registerStorage[state]);
+				accum += std::norm(registerStorage(i));
+				if (prob <= accum)
+				{
+					state = i;
+					break;
+				}
 			}
 
 			size_t measuredState = 0ULL;
 			size_t measuredStateMask = 0ULL;
-			if (prob <= accum)
+			if ((state & measuredQubitMask) != 0)
 			{
 				measuredState = 1ULL;
 				measuredStateMask = measuredQubitMask;
@@ -323,7 +328,6 @@ namespace QC {
 			const double norm = 1. / sqrt(accum);
 
 			// collapse
-
 			for (size_t state = 0; state < NrBasisStates; ++state)
 				registerStorage[state] *= ((state & measuredQubitMask) == measuredStateMask) ? norm : 0;
 
@@ -337,16 +341,24 @@ namespace QC {
 
 			const size_t measuredQubitMask = 1ULL << qubit;
 
-#pragma omp parallel for reduction(+:accum) num_threads(processor_count) schedule(static, OneQubitOmpLimit / divSchedule) 
-			for (long long state = measuredQubitMask; state < static_cast<long long int>(NrBasisStates); ++state)
+			size_t state = 0;
+			bool found = false;
+
+#pragma omp parallel for reduction(+:accum) num_threads(processor_count) schedule(static, OneQubitOmpLimit / divSchedule)
+			for (long long i = 0; i < NrBasisStates; ++i)
 			{
-				if (state & measuredQubitMask)
-					accum += std::norm(registerStorage[state]);
+				if (found) continue;
+				accum += std::norm(registerStorage(i));
+				if (prob <= accum)
+				{
+					state = i;
+					found = true;
+				}
 			}
-
+			
 			size_t measuredState = 0ULL;
 			size_t measuredStateMask = 0ULL;
-			if (prob <= accum)
+			if ((state & measuredQubitMask) != 0)
 			{
 				measuredState = 1ULL;
 				measuredStateMask = measuredQubitMask;
@@ -377,17 +389,18 @@ namespace QC {
 
 			const size_t measuredQubitMask = 1ULL << qubit;
 
-			for (size_t state = measuredQubitMask; state < NrBasisStates; ++state)
+			size_t state = 0;
+			for (size_t i = 0; i < NrBasisStates; ++i)
 			{
-				if (state & measuredQubitMask)
-					accum += std::norm(registerStorage[state]);
+				accum += std::norm(registerStorage(i));
+				if (prob <= accum)
+				{
+					state = i;
+					break;
+				}
 			}
 
-			size_t measuredState = 0;
-			if (prob <= accum)
-				measuredState = 1;
-
-			return measuredState;
+			return (state & measuredQubitMask) != 0 ? 1 : 0;
 		}
 
 		static inline size_t MeasureQubitNoCollapseOmp(size_t NrBasisStates, VectorClass& registerStorage, size_t qubit, const double prob)
@@ -397,18 +410,22 @@ namespace QC {
 
 			const size_t measuredQubitMask = 1ULL << qubit;
 
+			size_t state = 0;
+			bool found = false;
+
 #pragma omp parallel for reduction(+:accum) num_threads(processor_count) schedule(static, OneQubitOmpLimit / divSchedule)
-			for (long long state = measuredQubitMask; state < static_cast<long long int>(NrBasisStates); ++state)
+			for (long long i = 0; i < NrBasisStates; ++i)
 			{
-				if (state & measuredQubitMask)
-					accum += std::norm(registerStorage[state]);
+				if (found) continue;
+				accum += std::norm(registerStorage(i));
+				if (prob <= accum)
+				{
+					state = i;
+					found = true;
+				}
 			}
 
-			size_t measuredState = 0;
-			if (prob <= accum)
-				measuredState = 1;
-
-			return measuredState;
+			return (state & measuredQubitMask) != 0 ? 1 : 0;
 		}
 
 		static inline double GetQubitProbability(size_t NrBasisStates, const VectorClass& registerStorage, size_t qubit)
@@ -443,109 +460,64 @@ namespace QC {
 			return accum;
 		}
 
-		// TODO: the following are awful, there might be a better way to optimize this based on cache locality
-		// see the above that implement a single qubit measurement
 		static inline size_t Measure(size_t NrBasisStates, VectorClass& registerStorage, size_t firstQubit, size_t secondQubit, const double prob)
 		{
 			double accum = 0;
 
 			const size_t secondQubitp1 = secondQubit + 1;
-
 			const size_t firstPartMask = (1ULL << firstQubit) - 1;
 			const size_t measuredPartMask = (1ULL << secondQubitp1) - 1 - firstPartMask;
 			const size_t secondPartMask = NrBasisStates - 1 - measuredPartMask - firstPartMask;
 
-			const size_t secondPartMax = secondPartMask >> secondQubitp1;
-			const size_t maxMeasuredState = measuredPartMask >> firstQubit;
-
-			size_t measuredState = maxMeasuredState;
-
-			double norm = 1;
-			for (size_t state = 0; state <= maxMeasuredState; ++state)
+			size_t measuredState = 0;
+			for (size_t i = 0; i < NrBasisStates; ++i)
 			{
-				const size_t stateRegBits = state << firstQubit;
-				double stateProbability = 0;
-
-				for (size_t secondPartBits = 0; secondPartBits <= secondPartMax; ++secondPartBits)
-				{
-					const size_t secondPart = (secondPartBits << secondQubitp1) | stateRegBits;
-					for (size_t firstPartBits = 0; firstPartBits <= firstPartMask; ++firstPartBits)
-						stateProbability += std::norm(registerStorage[secondPart | firstPartBits]);
-				}
-
-				accum += stateProbability;
+				accum += std::norm(registerStorage(i));
 				if (prob <= accum)
 				{
-					measuredState = state;
-					norm = 1. / sqrt(stateProbability);
+					measuredState = i;
 					break;
 				}
 			}
+			measuredState &= measuredPartMask;
 
-			// collapse
-			for (size_t state = 0; state <= maxMeasuredState; ++state)
+			// find the norm
+			accum = 0;
+			for (size_t state = measuredState; state < NrBasisStates; ++state)
 			{
-				const size_t stateRegBits = state << firstQubit;
-
-				if (state == measuredState)
-				{
-					for (size_t secondPartBits = 0; secondPartBits <= secondPartMax; ++secondPartBits)
-					{
-						const size_t secondPart = (secondPartBits << secondQubitp1) | stateRegBits;
-						for (size_t firstPartBits = 0; firstPartBits <= firstPartMask; ++firstPartBits)
-							registerStorage[secondPart | firstPartBits] *= norm;
-					}
-				}
-				else
-				{
-					for (size_t secondPartBits = 0; secondPartBits <= secondPartMax; ++secondPartBits)
-					{
-						const size_t secondPart = (secondPartBits << secondQubitp1) | stateRegBits;
-						for (size_t firstPartBits = 0; firstPartBits <= firstPartMask; ++firstPartBits)
-							registerStorage[secondPart | firstPartBits] = 0;
-					}
-				}
+				if ((state & measuredPartMask) == measuredState)
+					accum += std::norm(registerStorage[state]);
 			}
+			const double norm = 1. / sqrt(accum);
 
+			// collapse
+			for (size_t state = 0; state < NrBasisStates; ++state)
+				registerStorage[state] *= ((state & measuredPartMask) == measuredState) ? norm : 0;
 
-			return measuredState;
+			return measuredState >> firstQubit;
 		}
 
 		static inline size_t MeasureNoCollapse(size_t NrBasisStates, VectorClass& registerStorage, size_t firstQubit, size_t secondQubit, const double prob)
 		{
 			double accum = 0;
 
 			const size_t secondQubitp1 = secondQubit + 1;
-
 			const size_t firstPartMask = (1ULL << firstQubit) - 1;
 			const size_t measuredPartMask = (1ULL << secondQubitp1) - 1 - firstPartMask;
 			const size_t secondPartMask = NrBasisStates - 1 - measuredPartMask - firstPartMask;
 
-			const size_t secondPartMax = secondPartMask >> secondQubitp1;
-			const size_t maxMeasuredState = measuredPartMask >> firstQubit;
-
-			size_t measuredState = maxMeasuredState;
-			for (size_t state = 0; state <= maxMeasuredState; ++state)
+			size_t measuredState = 0;
+			for (size_t i = 0; i < NrBasisStates; ++i)
 			{
-				const size_t stateRegBits = state << firstQubit;
-				double stateProbability = 0;
-
-				for (size_t secondPartBits = 0; secondPartBits <= secondPartMax; ++secondPartBits)
-				{
-					const size_t secondPart = (secondPartBits << secondQubitp1) | stateRegBits;
-					for (size_t firstPartBits = 0; firstPartBits <= firstPartMask; ++firstPartBits)
-						stateProbability += std::norm(registerStorage[secondPart | firstPartBits]);
-				}
-
-				accum += stateProbability;
+				accum += std::norm(registerStorage(i));
 				if (prob <= accum)
 				{
-					measuredState = state;
+					measuredState = i;
 					break;
 				}
 			}
 
-			return measuredState;
+			return (measuredState & measuredPartMask) >> firstQubit;
 		}
 
 		static int GetNumberOfThreads()

diff --git a/QCSim/Teleportation.h b/QCSim/Teleportation.h
@@ -144,7 +144,7 @@ namespace Teleportation
 		{
 			Teleport();
 
-			return AlgorithmClass::Measure(2, 2); // teleported qubit
+			return AlgorithmClass::Measure(2); // teleported qubit
 		}
 	};
-Original file line number
+Diff line change
@@ Expand Up / @@ -144,7 +144,7 @@ namespace Teleportation @@
     		{
     			Teleport();
-    			return AlgorithmClass::Measure(2, 2); // teleported qubit
+    			return AlgorithmClass::Measure(2); // teleported qubit
     		}
     	};
@@ Expand Down @@