added setDensityAmps

TysonRayJones · TysonRayJones · commit 8447c900b549 · 2021-10-14T02:19:51.000+11:00
diff --git a/QuEST/include/QuEST.h b/QuEST/include/QuEST.h
@@ -1494,7 +1494,8 @@ void initDebugState(Qureg qureg);
  */
 void initStateFromAmps(Qureg qureg, qreal* reals, qreal* imags);
 
-/** Overwrites a subset of the amplitudes in state-vector \p qureg, with those passed in \p reals and \p imags.
+/** Overwrites a contiguous subset of the amplitudes in state-vector \p qureg, 
+ * with those passed in \p reals and \p imags.
  *
  * Only amplitudes with indices in <b>[</b>\p startInd<b>,</b> \p startInd <b>+</b> \p numAmps<b>]</b> 
  * will be changed. The resulting \p qureg may not necessarily be in an L2 normalised state.
@@ -1522,6 +1523,7 @@ void initStateFromAmps(Qureg qureg, qreal* reals, qreal* imags);
  *
  *
  * @see
+ * - setDensityAmps()
  * - setWeightedQureg()
  * - initStateFromAmps()
  * - initBlankState()
@@ -1541,6 +1543,38 @@ void initStateFromAmps(Qureg qureg, qreal* reals, qreal* imags);
  */
 void setAmps(Qureg qureg, long long int startInd, qreal* reals, qreal* imags, long long int numAmps);
 
+/** Overwrites a contiguous subset of the amplitudes in density-matrix \p qureg, 
+ * with those passed in \p reals and \p imags, intrepreted column-wise.
+ *
+ * Only the first \p numAmp amplitudes starting from row-column index (\p startRow, \p startCol), and
+ * proceeding down the column (wrapping around between rows) will be modified. 
+ * The resulting \p qureg may not necessarily be a valid density matrix normalisation.
+ *
+ * In distributed mode, this function assumes the subset \p reals and \p imags exist
+ * (at least) on the node(s) containing the ultimately updated elements.\n
+ *
+ *
+ * @see
+ * - setAmps()
+ * - initStateFromAmps()
+ *
+ * @ingroup init
+ * @param[in,out] qureg the density-matrix to modify
+ * @param[in] startRow the row-index of the first amplitude in \p qureg to modify
+ * @param[in] startCol the column-index of the first amplitude in \p qureg to modify
+ * @param[in] reals array of the real components of the new amplitudes
+ * @param[in] imags array of the imaginary components of the new amplitudes
+ * @param[in] numAmps the length of each of the reals and imags arrays
+ * @throws invalidQuESTInputError()
+ * - if \p qureg is not a density matrix (i.e. is a state-vector)
+ * - if \p startRow is outside [0, `1 << qureg.numQubitsRepresented`]
+ * - if \p startCol is outside [0, `1 << qureg.numQubitsRepresented`]
+ * - if \p numAmps is outside [0, `qureg.numAmpsTotal`]
+ * - if \p numAmps is larger than the remaining number of amplitudes from (`startRow`, `startCol`), column-wise
+ * @author Tyson Jones
+ */
+void setDensityAmps(Qureg qureg, long long int startRow, long long int startCol, qreal* reals, qreal* imags, long long int numAmps);
+
 /** Overwrite the amplitudes of \p targetQureg with those from \p copyQureg. 
  * 
  * Registers must either both be state-vectors, or both be density matrices, and 
@@ -6621,7 +6655,7 @@ void applyQFT(Qureg qureg, int* qubits, int numQubits);
  * @ingroup operator
  * @param[in,out] qureg a state-vector or density matrix to modify
  * @param[in] qubit the qubit to which to apply the projector 
- * @param[in] the single-qubit outcome (`0` or `1`) to project \p qubit into
+ * @param[in] outcome the single-qubit outcome (`0` or `1`) to project \p qubit into
  * @throws invalidQuESTInputError()
  * - if \p qubit is outside [0, `qureg.numQubitsRepresented`)
  * - if \p outcome is not in {0,1}
diff --git a/QuEST/src/QuEST.c b/QuEST/src/QuEST.c
@@ -1027,9 +1027,12 @@ void setAmps(Qureg qureg, long long int startInd, qreal* reals, qreal* imags, lo
     qasm_recordComment(qureg, "Here, some amplitudes in the statevector were manually edited.");
 }
 
-void setDensityAmps(Qureg qureg, qreal* reals, qreal* imags) {
-    long long int numAmps = qureg.numAmpsTotal; 
-    statevec_setAmps(qureg, 0, reals, imags, numAmps);
+void setDensityAmps(Qureg qureg, long long int startRow, long long int startCol, qreal* reals, qreal* imags, long long int numAmps) {
+    validateDensityMatrQureg(qureg, __func__);
+    validateNumDensityAmps(qureg, startRow, startCol, numAmps, __func__);
+
+    long long int startInd = startRow + startCol*(1 << qureg.numQubitsRepresented);
+    statevec_setAmps(qureg, startInd, reals, imags, numAmps);
     
     qasm_recordComment(qureg, "Here, some amplitudes in the density matrix were manually edited.");
 }
diff --git a/QuEST/src/QuEST_validation.c b/QuEST/src/QuEST_validation.c
@@ -122,9 +122,9 @@ static const char* errorMessages[] = {
     [E_INVALID_STATE_INDEX] = "Invalid state index. Must be >=0 and <2^numQubits.",
     [E_INVALID_AMP_INDEX] = "Invalid amplitude index. Must be >=0 and <2^numQubits.",
     [E_INVALID_ELEM_INDEX] = "Invalid element index. Must be >=0 and <2^numQubits.",
-    [E_INVALID_NUM_AMPS] = "Invalid number of amplitudes. Must be >=0 and <=2^numQubits.",
+    [E_INVALID_NUM_AMPS] = "Invalid number of amplitudes. Must be >=0 and <=2^numQubits (or for density matrices, <=2^(2 numQubits)).",
     [E_INVALID_NUM_ELEMS] = "Invalid number of elements. Must be >=0 and <=2^numQubits.",
-    [E_INVALID_OFFSET_NUM_AMPS_QUREG] = "More amplitudes given than exist in the statevector from the given starting index.",
+    [E_INVALID_OFFSET_NUM_AMPS_QUREG] = "More amplitudes given than exist in the state from the given starting index.",
     [E_INVALID_OFFSET_NUM_ELEMS_DIAG] = "More elements given than exist in the diagonal operator from the given starting index.",
     [E_TARGET_IS_CONTROL] = "Control qubit cannot equal target qubit.",
     [E_TARGET_IN_CONTROLS] = "Control qubits cannot include target qubit.",
@@ -390,6 +390,15 @@ void validateNumAmps(Qureg qureg, long long int startInd, long long int numAmps,
     QuESTAssert(numAmps + startInd <= qureg.numAmpsTotal, E_INVALID_OFFSET_NUM_AMPS_QUREG, caller);
 }
 
+void validateNumDensityAmps(Qureg qureg, long long int startRow, long long int startCol, long long int numAmps, const char* caller) {
+    validateAmpIndex(qureg, startRow, caller);
+    validateAmpIndex(qureg, startCol, caller);
+    QuESTAssert(numAmps >= 0 && numAmps <= qureg.numAmpsTotal, E_INVALID_NUM_AMPS, caller);
+    
+    long long int startInd = startRow + startCol*(1 << qureg.numQubitsRepresented);
+    QuESTAssert(numAmps + startInd <= qureg.numAmpsTotal, E_INVALID_OFFSET_NUM_AMPS_QUREG, caller);
+}
+
 void validateNumElems(DiagonalOp op, long long int startInd, long long int numElems, const char* caller) {
     long long int indMax = 1LL << op.numQubits;
     QuESTAssert(startInd >= 0 && startInd < indMax, E_INVALID_ELEM_INDEX, caller);
diff --git a/QuEST/src/QuEST_validation.h b/QuEST/src/QuEST_validation.h
@@ -78,6 +78,8 @@ void validateSecondQuregStateVec(Qureg qureg2, const char *caller);
 
 void validateNumAmps(Qureg qureg, long long int startInd, long long int numAmps, const char* caller);
 
+void validateNumDensityAmps(Qureg qureg, long long int startRow, long long int startCol, long long int numAmps, const char* caller);
+
 void validateFileOpened(int opened, char* fn, const char* caller);
 
 void validateProb(qreal prob, const char* caller);
diff --git a/tests/test_state_initialisations.cpp b/tests/test_state_initialisations.cpp
@@ -444,6 +444,89 @@ TEST_CASE( "setAmps", "[state_initialisations]" ) {
 
 
 
+/** @sa setDensityAmps
+ * @ingroup unittest 
+ * @author Tyson Jones 
+ */
+TEST_CASE( "setDensityAmps", "[state_initialisations]" ) {
+    
+    Qureg matr = createDensityQureg(NUM_QUBITS, QUEST_ENV);
+    qreal reals[matr.numAmpsTotal];
+    qreal imags[matr.numAmpsTotal];
+    int maxInd = (1 << NUM_QUBITS);
+    
+    SECTION( "correctness" ) {
+        
+        SECTION( "density-matrix" ) {
+            
+            // all valid number of amplitudes and offsets
+            int startRow = GENERATE_COPY( range(0,maxInd) );
+            int startCol = GENERATE_COPY( range(0,maxInd) );
+            
+            int numPriorAmps = startRow + startCol*(1 << matr.numQubitsRepresented);
+            int maxNumAmps = matr.numAmpsTotal - numPriorAmps;
+            int numAmps = GENERATE_COPY( range(0,maxNumAmps) ); // upper-bound allows all amps specified
+            
+            // generate random amplitudes
+            for (int i=0; i<numAmps; i++) {
+                reals[i] = getRandomReal(-5,5);
+                imags[i] = getRandomReal(-5,5);
+            }
+            
+            // check both specified and un-specified amplitudes are correctly handled
+            initDebugState(matr);
+            QMatrix matrRef = toQMatrix(matr);
+            
+            setDensityAmps(matr, startRow, startCol, reals, imags, numAmps);
+            
+            int r=startRow;
+            int c=startCol;
+            for (int i=0; i<numAmps; i++) {
+                qcomp amp = reals[i] + imags[i] * (qcomp) 1i;
+                matrRef[r][c] = amp;
+                
+                r++;
+                if (r >= maxInd ) {
+                    r=0;
+                    c++;
+                }
+            }
+                
+            REQUIRE( areEqual(matr, matrRef) );
+        }
+    }
+    SECTION( "input validation" ) {
+        
+        SECTION( "start index" ) {
+            
+            int badInd = GENERATE_COPY( -1, maxInd );
+            int numAmps = 0;
+            REQUIRE_THROWS_WITH( setDensityAmps(matr, badInd, 0, reals, imags, numAmps), Contains("Invalid amplitude index") );
+            REQUIRE_THROWS_WITH( setDensityAmps(matr, 0, badInd, reals, imags, numAmps), Contains("Invalid amplitude index") );
+        }
+        
+        SECTION( "number of amplitudes" ) {
+            
+            // independent
+            int numAmps = GENERATE_COPY( -1, matr.numAmpsTotal+1 );
+            REQUIRE_THROWS_WITH( setDensityAmps(matr, 0, 0, reals, imags, numAmps), Contains("Invalid number of amplitudes") );
+
+            // invalid considering start-index
+            REQUIRE_THROWS_WITH( setDensityAmps(matr, maxInd-1, maxInd-1, reals, imags, 2), Contains("More amplitudes given than exist") );
+            REQUIRE_THROWS_WITH( setDensityAmps(matr, maxInd-1, maxInd-2, reals, imags, maxInd+2), Contains("More amplitudes given than exist") );
+        }
+        SECTION( "state-vector" ) {
+            
+            Qureg vec = createQureg(NUM_QUBITS, QUEST_ENV);
+            REQUIRE_THROWS_WITH( setDensityAmps(vec, 0, 0, reals, imags, 0), Contains("valid only for density matrices") );
+            destroyQureg(vec, QUEST_ENV);
+        }
+    }
+    destroyQureg(matr, QUEST_ENV);
+}
+
+
+
 /** @sa setWeightedQureg
  * @ingroup unittest 
  * @author Tyson Jones 

Original file line number	Diff line number	Diff line change
`@@ -1027,9 +1027,12 @@ void setAmps(Qureg qureg, long long int startInd, qreal* reals, qreal* imags, lo`
`1027`	`1027`	`qasm_recordComment(qureg, "Here, some amplitudes in the statevector were manually edited.");`
`1028`	`1028`	`}`
`1029`	`1029`
`1030`		`-void setDensityAmps(Qureg qureg, qreal* reals, qreal* imags) {`
`1031`		`- long long int numAmps = qureg.numAmpsTotal;`
`1032`		`- statevec_setAmps(qureg, 0, reals, imags, numAmps);`
	`1030`	`+void setDensityAmps(Qureg qureg, long long int startRow, long long int startCol, qreal* reals, qreal* imags, long long int numAmps) {`
	`1031`	`+ validateDensityMatrQureg(qureg, __func__);`
	`1032`	`+ validateNumDensityAmps(qureg, startRow, startCol, numAmps, __func__);`
	`1033`	`+`
	`1034`	`+ long long int startInd = startRow + startCol*(1 << qureg.numQubitsRepresented);`
	`1035`	`+ statevec_setAmps(qureg, startInd, reals, imags, numAmps);`
`1033`	`1036`
`1034`	`1037`	`qasm_recordComment(qureg, "Here, some amplitudes in the density matrix were manually edited.");`
`1035`	`1038`	`}`