diff --git a/tests/unit/trotterisation.cpp b/tests/unit/trotterisation.cpp index 35af81ba..883e6ed5 100644 --- a/tests/unit/trotterisation.cpp +++ b/tests/unit/trotterisation.cpp @@ -60,6 +60,47 @@ void TEST_ON_CACHED_QUREGS(quregCache quregs, auto& refFunc, auto& regularFunc, } } +void TEST_ON_CACHED_QUREGS(quregCache quregs, qvector& referenceResult, auto& testFunction, PauliStrSum& testHamiltonian) { + for (auto& [label, qureg]: quregs) { + + DYNAMIC_SECTION( label ) { + testFunction(qureg, testHamiltonian); + REQUIRE_AGREE(qureg, referenceResult); + } + + } + + return; +} + +void TEST_ON_CACHED_QUREGS(quregCache quregs, qmatrix& referenceResult, auto& testFunction, PauliStrSum& testHamiltonian) { + for (auto& [label, qureg]: quregs) { + + DYNAMIC_SECTION( label ) { + testFunction(qureg, testHamiltonian); + REQUIRE_AGREE(qureg, referenceResult); + } + + } + + return; +} + +void TEST_OBSERVABLES_ON_QUREGS(quregCache quregs, qvector& referenceResult, auto& testFunction, PauliStrSum testHamiltonian, PauliStrSum testObservable) { + for (auto& [label, qureg]: quregs) { + + DYNAMIC_SECTION( label ) { + qvector testResult = testFunction(qureg, testHamiltonian, testObservable); + REQUIRE_AGREE(calcTotalProb(qureg), 1.0); + REQUIRE_AGREE(testResult, referenceResult); + } + + } + + return; +} + + /* * Prepare a Hamiltonian H under which dynamical * evolution will be simulated via Trotterisation @@ -220,29 +261,7 @@ TEST_CASE( "randomisedTrotter", TEST_CATEGORY ) { */ TEST_CASE( "applyTrotterizedUnitaryTimeEvolution", TEST_CATEGORY ) { - // BEWARE: this test creates a new Qureg below which will have - // deployments chosen by the auto-deployer; it is ergo unpredictable - // whether it will be multithreaded, GPU-accelerated or distributed. - // This test is ergo checking only a single, unspecified deployment, - // unlike other tests which check all deployments. This is tolerable - // since (non-randomised) Trotterisation is merely invoking routines - // (Pauli gadgets) already independently tested across deployments - SECTION( LABEL_CORRECTNESS ) { - - int numQubits = 20; - Qureg qureg = createQureg(numQubits); - initPlusState(qureg); - bool permutePaulis = false; - - PauliStrSum hamil = createHeisenbergHamiltonian(numQubits); - PauliStrSum observ = createAlternatingPauliObservable(numQubits); - - qreal dt = 0.1; - int order = 4; - int reps = 5; - int steps = 10; - // nudge the epsilon used by internal validation functions up a bit // as the time evolution operation plays badly with single precision // Defaults for validation epsilon are: @@ -252,57 +271,68 @@ TEST_CASE( "applyTrotterizedUnitaryTimeEvolution", TEST_CATEGORY ) { qreal initialValidationEps = getValidationEpsilon(); setValidationEpsilon(2 * initialValidationEps); - /* - * Tolerance for floating-point comparisons - * Note that the underlying numerics are sensitive to the float - * precision AND to the number of threads. As such we set quite - * large epsilon values to account for the worst-case scenario which - * is single precision, single thread. The baseline for these results - * is double precision, multiple threads. - * - * Values (assuming default initialValidationEps) are: - * Single precision: - * obsEps = 0.03 - * normEps = 0.001 - * - * Double precision: - * obsEps = 3E-9 - * normEps = 1E-10 - * - * Quad precision: - * obsEps = 3E-10 - * normEps = 1E-11 - */ - qreal obsEps = 3E3 * initialValidationEps; - qreal normEps = 100 * initialValidationEps; - - vector refObservables = { - 19.26827777028073, - 20.34277275871839, - 21.21120737889526, - 21.86585902741717, - 22.30371711358924, - 22.52644660547882, - 22.54015748825067, - 22.35499202583118, - 21.9845541501027, - 21.44521638719462 + const int NUM_QUBITS = 8; + qreal dt = 0.1; + int order = 4; + int reps = 5; + const int STEPS = 20; + bool permutePaulis = GENERATE(true, false); + + auto unitaryTimeEvoFunc = + [dt, order, reps, STEPS, permutePaulis](Qureg qureg, PauliStrSum& hamil, PauliStrSum& observable) + -> qvector { + qvector observations = getZeroVector(STEPS); + initPlusState(qureg); + + for (int i = 0; i < STEPS; i++) { + applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps, permutePaulis); + observations.at(i) = calcExpecPauliStrSum(qureg, observable); + } + + return observations; }; - for (int i = 0; i < steps; i++) { - applyTrotterizedUnitaryTimeEvolution(qureg, hamil, dt, order, reps, permutePaulis); - qreal expec = calcExpecPauliStrSum(qureg, observ); - - REQUIRE_THAT( expec, WithinAbs(refObservables[i], obsEps) ); - } + PauliStrSum hamil = createHeisenbergHamiltonian(NUM_QUBITS); + PauliStrSum observ = createAlternatingPauliObservable(NUM_QUBITS); - // Verify state remains normalized - REQUIRE_THAT( calcTotalProb(qureg), WithinAbs(1.0, normEps) ); + qvector refObservables = { + 8.521995598825049, + 8.963711845322115, + 9.32005226684505, + 9.587768088649602, + 9.765522600223822, + 9.85387668440598, + 9.855195944206464, + 9.773484879367675, + 9.614158409472378, + 9.383765238225045, + 9.089680663909942, + 8.739788123639109, + 8.342168826039893, + 7.904817272753528, + 7.435397472039873, + 6.94105054616863, + 6.428259679798389, + 5.90277345392904, + 5.369584051930907, + 4.832953030744839 + }; + + SECTION("Time Evolve Statevectors") { + quregCache eightQubitSVCache = createFixedSizeCachedStatevecsOrDensmatrs(NUM_QUBITS, false); + TEST_OBSERVABLES_ON_QUREGS(eightQubitSVCache, refObservables, unitaryTimeEvoFunc, hamil, observ); + destroyCache(eightQubitSVCache); + } + + SECTION("Time Evolve Density Matrices") { + quregCache eightQubitDMCache = createFixedSizeCachedStatevecsOrDensmatrs(NUM_QUBITS, true); + TEST_OBSERVABLES_ON_QUREGS(eightQubitDMCache, refObservables, unitaryTimeEvoFunc, hamil, observ); + destroyCache(eightQubitDMCache); + } // Restore validation epsilon setValidationEpsilon(initialValidationEps); - destroyQureg(qureg); destroyPauliStrSum(hamil); destroyPauliStrSum(observ); } @@ -404,146 +434,122 @@ TEST_CASE( "applyTrotterizedUnitaryTimeEvolution", TEST_CATEGORY ) { TEST_CASE( "applyTrotterizedImaginaryTimeEvolution", TEST_CATEGORY ) { - - // BEWARE: this test creates a new Qureg below which will have - // deployments chosen by the auto-deployer; it is ergo unpredictable - // whether it will be multithreaded, GPU-accelerated or distributed. - // This test is ergo checking only a single, unspecified deployment, - // unlike other tests which check all deployments. This is tolerable - // since (non-randomised) Trotterisation is merely invoking routines - // (Pauli gadgets) already independently tested across deployments + int numQubits = getNumCachedQubits(); + auto statevecQuregs = getCachedStatevecs(); + auto densmatrQuregs = getCachedDensmatrs(); SECTION( LABEL_CORRECTNESS ) { - int numQubits = 16; qreal tau = 0.1; int order = 6; int reps = 5; int steps = 10; - bool permutePaulis = false; - - // Tolerance for ground state amplitude - qreal eps = 1E-2; + bool permutePaulis = GENERATE(true, false); + + auto driveToGroundFunc = [steps, tau, order, reps, permutePaulis](Qureg qureg, PauliStrSum& hamil) { + initPlusState(qureg); + for (int i = 0; i < steps; ++i) { + applyTrotterizedImaginaryTimeEvolution(qureg, hamil, tau, order, reps, permutePaulis); + setQuregToRenormalized(qureg); + } + }; + + +#if FLOAT_PRECISION == 4 + /* + * The numerical exponent is sufficiently inaccurate to breach the default + * tolerances at quad precision, so we apply the following kludge to prevent irritating test failures. + * The real lessons from these tests are: + * - Don't do time-evolution at single precision. + * - Don't do time-evolution in serial. + */ + + qreal initialEps = getTestAbsoluteEpsilon(); + setTestAbsoluteEpsilon(30 * initialEps); +#endif + + // Ground state: all qubits align down (driven by strong magnetic field) + SECTION("Spin Down Field") { - Qureg qureg = createQureg(numQubits); - initPlusState(qureg); - PauliStrSum ising = createIsingHamiltonian(numQubits, 10.0, 0.0, 0.0); - - for (int i = 0; i < steps; ++i) { - applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis); - setQuregToRenormalized(qureg); - } - - qcomp amp = getQuregAmp(qureg, 0); - qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag(); - - REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) ); - - for (long long i = 1; i < (1LL << numQubits); i++) { - qcomp other_amp = getQuregAmp(qureg, i); - qreal other_mag = other_amp.real() * other_amp.real() + - other_amp.imag() * other_amp.imag(); - REQUIRE( other_mag < eps ); - } - - destroyQureg(qureg); + + qvector statevecRef = getZeroVector(getPow2(numQubits)); + statevecRef.at(0) = 1; + + qmatrix densmatrRef = getZeroMatrix(getPow2(numQubits)); + densmatrRef[0][0] = 1; + + TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, driveToGroundFunc, ising); + TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, driveToGroundFunc, ising); + destroyPauliStrSum(ising); } // Ground state: all qubits align up (driven by opposite magnetic field) - { - Qureg qureg = createQureg(numQubits); - initPlusState(qureg); - + SECTION("Spin Up Field") + { PauliStrSum ising = createIsingHamiltonian(numQubits, -10.0, 0.0, 0.0); - - for (int i = 0; i < steps; ++i) { - applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis); - setQuregToRenormalized(qureg); - } - - long long last_state = (1LL << numQubits) - 1; - qcomp amp = getQuregAmp(qureg, last_state); - qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag(); - - REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) ); - - for (long long i = 0; i < (1LL << numQubits); i++) { - if (i == last_state) continue; - qcomp other_amp = getQuregAmp(qureg, i); - qreal other_mag = other_amp.real() * other_amp.real() + - other_amp.imag() * other_amp.imag(); - REQUIRE( other_mag < eps ); - } - - destroyQureg(qureg); + + qindex namps = getPow2(numQubits); + + qvector statevecRef = getZeroVector(namps); + statevecRef.at(namps - 1) = 1; + + qmatrix densmatrRef = getZeroMatrix(namps); + densmatrRef[namps-1][namps-1] = 1; + + TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, driveToGroundFunc, ising); + TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, driveToGroundFunc, ising); + destroyPauliStrSum(ising); } // Ground state: all qubits align down (driven by ferromagnetic interactions and bias) + SECTION("Ferromagnetic Interaction") { - Qureg qureg = createQureg(numQubits); - initPlusState(qureg); - PauliStrSum ising = createIsingHamiltonian(numQubits, 0.0, 10.0, 10.0); - - for (int i = 0; i < steps; ++i) { - applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis); - setQuregToRenormalized(qureg); - } - - qcomp amp = getQuregAmp(qureg, 0); - qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag(); - - REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) ); - - for (long long i = 1; i < (1LL << numQubits); i++) { - qcomp other_amp = getQuregAmp(qureg, i); - qreal other_mag = other_amp.real() * other_amp.real() + - other_amp.imag() * other_amp.imag(); - REQUIRE( other_mag < eps ); - } - - destroyQureg(qureg); + + qvector statevecRef = getZeroVector(getPow2(numQubits)); + statevecRef.at(0) = 1; + + qmatrix densmatrRef = getZeroMatrix(getPow2(numQubits)); + densmatrRef[0][0] = 1; + + TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, driveToGroundFunc, ising); + TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, driveToGroundFunc, ising); + destroyPauliStrSum(ising); } // Ground state: alternating pattern (driven by antiferromagnetic interactions) + SECTION("Antiferromagnetic Interaction") { - Qureg qureg = createQureg(numQubits); - initPlusState(qureg); - PauliStrSum ising = createIsingHamiltonian(numQubits, 0.0, -10.0, 10.0); - for (int i = 0; i < steps; ++i) { - applyTrotterizedImaginaryTimeEvolution(qureg, ising, tau, order, reps, permutePaulis); - setQuregToRenormalized(qureg); - } - + // This should correctly pick out the non-zero amplitude + // Qubit 0 is always 0 thanks to asymmetric bias unsigned long long idx = 0; for (int i = 0; i < numQubits / 2; ++i) { idx += (1ULL << (2*i + 1)); } - - qcomp amp = getQuregAmp(qureg, idx); - qreal amp_mag = amp.real() * amp.real() + amp.imag() * amp.imag(); - - REQUIRE_THAT( amp_mag, WithinAbs(1.0, eps) ); - - for (long long i = 0; i < (1LL << numQubits); i++) { - if (i == idx) continue; - qcomp other_amp = getQuregAmp(qureg, i); - qreal other_mag = other_amp.real() * other_amp.real() + - other_amp.imag() * other_amp.imag(); - REQUIRE( other_mag < eps ); - } - - destroyQureg(qureg); + + qvector statevecRef = getZeroVector(getPow2(numQubits)); + statevecRef.at(idx) = 1; + + qmatrix densmatrRef = getZeroMatrix(getPow2(numQubits)); + densmatrRef[idx][idx] = 1; + + TEST_ON_CACHED_QUREGS(statevecQuregs, statevecRef, driveToGroundFunc, ising); + TEST_ON_CACHED_QUREGS(densmatrQuregs, densmatrRef, driveToGroundFunc, ising); + destroyPauliStrSum(ising); } + +#if FLOAT_PRECISION == 4 + setTestEpsilon(initialEps); +#endif } SECTION( LABEL_VALIDATION ) { diff --git a/tests/utils/cache.cpp b/tests/utils/cache.cpp index 0cb8a603..186d8365 100644 --- a/tests/utils/cache.cpp +++ b/tests/utils/cache.cpp @@ -91,17 +91,21 @@ deployInfo getSupportedDeployments() { * manage cached quregs */ -quregCache createCachedStatevecsOrDensmatrs(bool isDensMatr) { +quregCache createFixedSizeCachedStatevecsOrDensmatrs(const int NUM_QUBITS, const bool IS_DENSITY_MATRIX) { quregCache out; // only add supported-deployment quregs to the cache for (auto [label, mpi, gpu, omp] : getSupportedDeployments()) - out[label] = createCustomQureg(getNumCachedQubits(), isDensMatr, mpi, gpu, omp); + out[label] = createCustomQureg(NUM_QUBITS, IS_DENSITY_MATRIX, mpi, gpu, omp); return out; } +quregCache createCachedStatevecsOrDensmatrs(const bool IS_DENSITY_MATRIX) { + return createFixedSizeCachedStatevecsOrDensmatrs(getNumCachedQubits(), IS_DENSITY_MATRIX); +} + void createCachedQuregs() { // must not be called twice nor pre-creation @@ -116,6 +120,12 @@ void createCachedQuregs() { densmatrs2 = createCachedStatevecsOrDensmatrs(true); } +void destroyCache(quregCache& cache) { + for (auto& [label, qureg]: cache) + destroyQureg(qureg); + return; +} + void destroyCachedQuregs() { // must not be called twice nor pre-creation @@ -124,13 +134,12 @@ void destroyCachedQuregs() { DEMAND( ! densmatrs1.empty() ); DEMAND( ! densmatrs2.empty() ); - auto caches = { + std::vector caches = { statevecs1, statevecs2, densmatrs1, densmatrs2}; for (auto& cache : caches) - for (auto& [label, qureg]: cache) - destroyQureg(qureg); + destroyCache(cache); statevecs1.clear(); statevecs2.clear(); diff --git a/tests/utils/cache.hpp b/tests/utils/cache.hpp index 87d8382f..36b1971d 100644 --- a/tests/utils/cache.hpp +++ b/tests/utils/cache.hpp @@ -29,11 +29,15 @@ using deployInfo = std::vector>; int getNumCachedQubits(); deployInfo getSupportedDeployments(); +quregCache createFixedSizeCachedStatevecsOrDensmatrs(const int NUM_QUBITS, const bool IS_DENSITY_MATRIX); +quregCache createCachedStatevecsOrDensmatrs(const bool IS_DENSITY_MATRIX); + void createCachedFullStateDiagMatrs(); void destroyCachedFullStateDiagMatrs(); matrixCache getCachedFullStateDiagMatrs(); void createCachedQuregs(); +void destroyCache(quregCache& cache); void destroyCachedQuregs(); quregCache getCachedStatevecs(); quregCache getCachedDensmatrs(); diff --git a/tests/utils/compare.cpp b/tests/utils/compare.cpp index c3558543..3af15d60 100644 --- a/tests/utils/compare.cpp +++ b/tests/utils/compare.cpp @@ -34,25 +34,42 @@ using namespace Catch::Matchers; #if FLOAT_PRECISION == 1 - const qreal ABSOLUTE_EPSILON = 1E-2; - const qreal RELATIVE_EPSILON = 1E-2; + qreal absolute_epsilon = 1E-2; + qreal relative_epsilon = 1E-2; #elif FLOAT_PRECISION == 2 - const qreal ABSOLUTE_EPSILON = 1E-8; - const qreal RELATIVE_EPSILON = 1E-8; + qreal absolute_epsilon = 1E-8; + qreal relative_epsilon = 1E-8; #elif FLOAT_PRECISION == 4 - const qreal ABSOLUTE_EPSILON = 1E-10; - const qreal RELATIVE_EPSILON = 1E-10; + qreal absolute_epsilon = 1E-10; + qreal relative_epsilon = 1E-10; #endif qreal getTestAbsoluteEpsilon() { - return ABSOLUTE_EPSILON; + return absolute_epsilon; } qreal getTestRelativeEpsilon() { - return RELATIVE_EPSILON; + return relative_epsilon; +} + + +void setTestAbsoluteEpsilon(const qreal TEST_ABS_EPS) { + absolute_epsilon = TEST_ABS_EPS; + return; +} + +void setTestRelativeEpsilon(const qreal TEST_REL_EPS) { + relative_epsilon = TEST_REL_EPS; + return; +} + +void setTestEpsilon(const qreal TEST_EPS) { + setTestAbsoluteEpsilon(TEST_EPS); + setTestRelativeEpsilon(TEST_EPS); + return; } @@ -72,10 +89,10 @@ bool doScalarsAgree(qcomp a, qcomp b) { // permit absolute OR relative agreement - if (getAbsDif(a, b) <= ABSOLUTE_EPSILON) + if (getAbsDif(a, b) <= absolute_epsilon) return true; - return (getRelDif(a, b) <= RELATIVE_EPSILON); + return (getRelDif(a, b) <= relative_epsilon); } bool doMatricesAgree(qmatrix a, qmatrix b) { @@ -107,8 +124,8 @@ void REPORT_AMP_AND_FAIL( size_t index, qcomp amplitude, qcomp reference ) { qreal relative_difference = getRelDif(amplitude, reference); CAPTURE( index, amplitude, reference, - absolute_difference, ABSOLUTE_EPSILON, - relative_difference, RELATIVE_EPSILON + absolute_difference, absolute_epsilon, + relative_difference, relative_epsilon ); FAIL( ); } @@ -154,8 +171,8 @@ void REPORT_SCALAR_AND_FAIL( qcomp scalar, qcomp reference ) { qreal relative_difference = getRelDif(scalar, reference); CAPTURE( scalar, reference, - absolute_difference, ABSOLUTE_EPSILON, - relative_difference, RELATIVE_EPSILON + absolute_difference, absolute_epsilon, + relative_difference, relative_epsilon ); FAIL( ); } @@ -167,8 +184,8 @@ void REPORT_SCALAR_AND_FAIL( qreal scalar, qreal reference ) { qreal relative_difference = getRelDif(qcomp(scalar,0), qcomp(reference,0)); CAPTURE( scalar, reference, - absolute_difference, ABSOLUTE_EPSILON, - relative_difference, RELATIVE_EPSILON + absolute_difference, absolute_epsilon, + relative_difference, relative_epsilon ); FAIL( ); } @@ -230,8 +247,8 @@ void REPORT_ELEM_AND_FAIL( size_t row, size_t col, qcomp elem, qcomp reference ) qreal relative_difference = getRelDif(elem, reference); CAPTURE( row, col, elem, reference, - absolute_difference, ABSOLUTE_EPSILON, - relative_difference, RELATIVE_EPSILON + absolute_difference, absolute_epsilon, + relative_difference, relative_epsilon ); FAIL( ); } diff --git a/tests/utils/compare.hpp b/tests/utils/compare.hpp index 028b7202..96d11f10 100644 --- a/tests/utils/compare.hpp +++ b/tests/utils/compare.hpp @@ -24,6 +24,10 @@ using std::vector; qreal getTestAbsoluteEpsilon(); qreal getTestRelativeEpsilon(); +void setTestAbsoluteEpsilon(const qreal TEST_ABS_EPS); +void setTestRelativeEpsilon(const qreal TEST_REL_EPS); +void setTestEpsilon(const qreal TEST_EPS); + bool doScalarsAgree(qcomp a, qcomp b); bool doMatricesAgree(qmatrix a, qmatrix b);