Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

🧹 Noise-Aware Simulator Cleanup #491

Merged
merged 31 commits into from
Feb 5, 2024
Merged

🧹 Noise-Aware Simulator Cleanup #491

Show file tree
Hide file tree
Changes from 22 commits
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
d53b81f
adding support for non-unitary functions for density matrix sim
Sep 8, 2023
df36884
Merge remote-tracking branch 'origin/main' into nw_sim_features
Oct 11, 2023
098e7b8
removing unused parameter.
Oct 16, 2023
c229337
cleanup
Oct 16, 2023
10398a6
adding support for non unitary operations in the density matrix simul…
Nov 30, 2023
3fc57d6
Merge remote-tracking branch 'origin/main' into nw_sim_features
Nov 30, 2023
7aec8a7
fix stuff after merge
Nov 30, 2023
6727fc6
some cleanup
Dec 1, 2023
7298bd1
Removing the last traces of sequentially apply noise
Dec 1, 2023
3de7e63
update tests
Dec 1, 2023
68f5ee9
update more tests
Dec 1, 2023
8e1c566
add a new test
Dec 1, 2023
88a3d13
🎨 pre-commit fixes
pre-commit-ci[bot] Dec 1, 2023
09668f4
addressing CI reports
Dec 1, 2023
282b785
removing unused variable
Dec 1, 2023
3abf2ab
Merge branch 'main' into nw_sim_features
Jan 12, 2024
3149d3b
Merge branch 'main' into nw_sim_features
Jan 12, 2024
53a393e
🎨 pre-commit fixes
pre-commit-ci[bot] Jan 19, 2024
2f18ecb
remove unused code
Jan 19, 2024
b03a1b6
adding new test
Jan 19, 2024
56658e2
🎨 pre-commit fixes
pre-commit-ci[bot] Jan 19, 2024
b07dae0
Merge branch 'main' into nw_sim_features
Jan 31, 2024
09fe61f
addressing comments
Jan 31, 2024
0f465d7
🎨 pre-commit fixes
pre-commit-ci[bot] Jan 31, 2024
5c5fce1
addressing comments
Feb 5, 2024
b24ce4c
merging changes from remote
Feb 5, 2024
eec502f
🎨 pre-commit fixes
pre-commit-ci[bot] Feb 5, 2024
9dc1a40
removing unused code
Feb 5, 2024
2bcc315
🎨 pre-commit fixes
pre-commit-ci[bot] Feb 5, 2024
b3527be
Merge branch 'main' into nw_sim_features
Feb 5, 2024
c71e92e
Merge branch 'main' into nw_sim_features
33Gjl1Xe Feb 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
192 changes: 12 additions & 180 deletions include/mqt-core/dd/NoiseFunctionality.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -254,16 +254,13 @@ template <class Config> class DeterministicNoiseFunctionality {
double noiseProbabilityMultiQubit,
double ampDampProbSingleQubit,
double ampDampProbMultiQubit,
std::vector<NoiseOperations> effects,
bool useDensityMatType, bool seqApplyNoise)
std::vector<NoiseOperations> effects)
: package(dd), nQubits(nq),
noiseProbSingleQubit(noiseProbabilitySingleQubit),
noiseProbMultiQubit(noiseProbabilityMultiQubit),
ampDampingProbSingleQubit(ampDampProbSingleQubit),
ampDampingProbMultiQubit(ampDampProbMultiQubit),
noiseEffects(std::move(effects)),
useDensityMatrixType(useDensityMatType),
sequentiallyApplyNoise(seqApplyNoise) {}
noiseEffects(std::move(effects)) {}

protected:
// NOLINTNEXTLINE(cppcoreguidelines-avoid-const-or-ref-data-members)
Expand All @@ -276,43 +273,23 @@ template <class Config> class DeterministicNoiseFunctionality {
double ampDampingProbMultiQubit;

std::vector<NoiseOperations> noiseEffects;
bool useDensityMatrixType;
bool sequentiallyApplyNoise;

inline static const std::map<NoiseOperations, std::size_t>
SEQUENTIAL_NOISE_MAP = {
{Identity, 1}, // Identity Noise
{PhaseFlip, 2}, // Phase-flip
{AmplitudeDamping, 2}, // Amplitude Damping
{Depolarization, 4}, // Depolarisation
};

[[nodiscard]] std::size_t getNumberOfQubits() const { return nQubits; }

public:
void applyNoiseEffects(dEdge& originalEdge,
const std::unique_ptr<qc::Operation>& qcOperation) {
auto usedQubits = qcOperation->getUsedQubits();
if (sequentiallyApplyNoise) {
applyDetNoiseSequential(originalEdge, usedQubits);
} else {
dCachedEdge nodeAfterNoise = {};
if (useDensityMatrixType) {
dEdge::applyDmChangesToEdge(originalEdge);
nodeAfterNoise = applyNoiseEffects(originalEdge, usedQubits, false);
dEdge::revertDmChangesToEdge(originalEdge);
} else {
nodeAfterNoise = applyNoiseEffects(originalEdge, usedQubits, true);
}
auto r = dEdge{nodeAfterNoise.p, package->cn.lookup(nodeAfterNoise.w)};
package->incRef(r);
dEdge::alignDensityEdge(originalEdge);
package->decRef(originalEdge);
originalEdge = r;
if (useDensityMatrixType) {
dEdge::setDensityMatrixTrue(originalEdge);
}
}
dCachedEdge nodeAfterNoise = {};
dEdge::applyDmChangesToEdge(originalEdge);
nodeAfterNoise = applyNoiseEffects(originalEdge, usedQubits, false);
dEdge::revertDmChangesToEdge(originalEdge);
auto r = dEdge{nodeAfterNoise.p, package->cn.lookup(nodeAfterNoise.w)};
package->incRef(r);
dEdge::alignDensityEdge(originalEdge);
package->decRef(originalEdge);
originalEdge = r;
dEdge::setDensityMatrixTrue(originalEdge);
}

private:
Expand Down Expand Up @@ -486,151 +463,6 @@ template <class Config> class DeterministicNoiseFunctionality {
e[3] = package->add2(helperEdge[0], helperEdge[1], var);
}
}

void applyDetNoiseSequential(dEdge& originalEdge,
const std::set<qc::Qubit>& targets) {
std::array<mEdge, NrEdges::value> idleOperation{};

// Iterate over qubits and check if the qubit had been used
for (const auto targetQubit : targets) {
for (auto const& type : noiseEffects) {
generateGate(idleOperation, type, targetQubit,
getNoiseProbability(type, targets));
std::optional<dEdge> tmp{};
// Apply all noise matrices of the current noise effect
for (std::size_t m = 0; m < SEQUENTIAL_NOISE_MAP.find(type)->second;
m++) {
auto tmp0 = package->conjugateTranspose(idleOperation.at(m));
auto tmp1 = package->multiply(originalEdge,
densityFromMatrixEdge(tmp0), 0, false);
auto tmp2 =
package->multiply(densityFromMatrixEdge(idleOperation.at(m)),
tmp1, 0, useDensityMatrixType);
if (!tmp.has_value()) {
tmp = tmp2;
} else {
tmp = package->add(tmp2, *tmp);
}
}
assert(tmp.has_value());
auto& tmpEdge = *tmp;
package->incRef(tmpEdge);
dEdge::alignDensityEdge(originalEdge);
package->decRef(originalEdge);
originalEdge = tmpEdge;
if (useDensityMatrixType) {
dEdge::setDensityMatrixTrue(originalEdge);
}
}
}
}

void generateDepolarizationGate(
std::array<mEdge, NrEdges::value>& pointerForMatrices,
const qc::Qubit target, const double probability) {
std::array<GateMatrix, NrEdges::value> idleNoiseGate{};

// (1 0)
// sqrt(1- ((3p)/4))*(0 1)
idleNoiseGate[0][0] = idleNoiseGate[0][3] =
std::sqrt(1 - ((3 * probability) / 4));
idleNoiseGate[0][1] = idleNoiseGate[0][2] = 0;

pointerForMatrices[0] =
package->makeGateDD(idleNoiseGate[0], getNumberOfQubits(), target);

// (0 1)
// sqrt(probability/4))*(1 0)
idleNoiseGate[1][1] = idleNoiseGate[1][2] = std::sqrt(probability / 4);
idleNoiseGate[1][0] = idleNoiseGate[1][3] = 0;

pointerForMatrices[1] =
package->makeGateDD(idleNoiseGate[1], getNumberOfQubits(), target);

// (1 0)
// sqrt(probability/4))*(0 -1)
idleNoiseGate[2][0] = std::sqrt(probability / 4);
idleNoiseGate[2][3] = -std::sqrt(probability / 4);
idleNoiseGate[2][1] = idleNoiseGate[2][2] = 0;

pointerForMatrices[3] =
package->makeGateDD(idleNoiseGate[2], getNumberOfQubits(), target);

// (0 -i)
// sqrt(probability/4))*(i 0)
idleNoiseGate[3][2] = {0, std::sqrt(probability / 4)};
idleNoiseGate[3][1] = {0, -std::sqrt(probability / 4)};
idleNoiseGate[3][0] = idleNoiseGate[3][3] = 0;

pointerForMatrices[2] =
package->makeGateDD(idleNoiseGate[3], getNumberOfQubits(), target);
}

void generateGate(std::array<mEdge, NrEdges::value>& pointerForMatrices,
const NoiseOperations noiseType, const qc::Qubit target,
const double probability) {
std::array<GateMatrix, NrEdges::value> idleNoiseGate{};
switch (noiseType) {
// identity noise (for testing)
// (1 0)
// (0 1),
case Identity: {
pointerForMatrices[0] =
package->makeGateDD(I_MAT, getNumberOfQubits(), target);

break;
}
// phase flip
// (1 0) (1 0)
// e0= sqrt(1-probability)*(0 1), e1= sqrt(probability)*(0 -1)
case PhaseFlip: {
idleNoiseGate[0][0] = idleNoiseGate[0][3] = std::sqrt(1 - probability);
idleNoiseGate[0][1] = idleNoiseGate[0][2] = 0;
idleNoiseGate[1][0] = std::sqrt(probability);
idleNoiseGate[1][3] = -std::sqrt(probability);
idleNoiseGate[1][1] = idleNoiseGate[1][2] = 0;

pointerForMatrices[0] =
package->makeGateDD(idleNoiseGate[0], getNumberOfQubits(), target);
pointerForMatrices[1] =
package->makeGateDD(idleNoiseGate[1], getNumberOfQubits(), target);

break;
}
// amplitude damping
// (1 0) (0 sqrt(probability))
// e0= (0 sqrt(1-probability), e1= (0 0)
case AmplitudeDamping: {
idleNoiseGate[0][0] = 1;
idleNoiseGate[0][1] = idleNoiseGate[0][2] = 0;
idleNoiseGate[0][3] = std::sqrt(1 - probability);

idleNoiseGate[1][0] = idleNoiseGate[1][3] = idleNoiseGate[1][2] = 0;
idleNoiseGate[1][1] = std::sqrt(probability);

pointerForMatrices[0] =
package->makeGateDD(idleNoiseGate[0], getNumberOfQubits(), target);
pointerForMatrices[1] =
package->makeGateDD(idleNoiseGate[1], getNumberOfQubits(), target);
break;
}
// depolarization
case Depolarization:
generateDepolarizationGate(pointerForMatrices, target, probability);
break;
default:
throw std::runtime_error("Unknown noise effect received.");
}
}

double getNoiseProbability(const NoiseOperations type,
const std::set<qc::Qubit>& targets) {
if (type == AmplitudeDamping) {
return (targets.size() == 1) ? ampDampingProbSingleQubit
: ampDampingProbMultiQubit;
}
return (targets.size() == 1) ? noiseProbSingleQubit : noiseProbMultiQubit;
}
};

} // namespace dd
103 changes: 82 additions & 21 deletions include/mqt-core/dd/Package.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1165,6 +1165,71 @@ template <class Config> class Package {
return '1';
}

mEdge buildMeasOp(const Qubit index, const char measureFor = '0') {
mEdge measOp = {};
mEdge f = {};

const auto identityMatrix =
std::array{mEdge::one(), mEdge::zero(), mEdge::zero(), mEdge::one()};

const auto measureMatrix = measureFor == '0'
? std::array{mEdge::one(), mEdge::zero(),
mEdge::zero(), mEdge::zero()}
: std::array{mEdge::zero(), mEdge::zero(),
mEdge::zero(), mEdge::one()};

if (index == 0) {
measOp = makeDDNode(static_cast<dd::Qubit>(0), measureMatrix);
} else {
measOp = makeDDNode(static_cast<dd::Qubit>(0), identityMatrix);
}

for (dd::Qubit p = 1; p < static_cast<dd::Qubit>(qubits()); p++) {
if (p == index) {
f = makeDDNode(static_cast<dd::Qubit>(0), measureMatrix);
} else {
f = makeDDNode(static_cast<dd::Qubit>(0), identityMatrix);
}
measOp = kronecker(f, measOp);
}
return measOp;
}

std::pair<dEdge, char> measureOneCollapsing(dEdge& e, const Qubit index,
std::mt19937_64& mt) {
char measuredResult = '0';

auto operation = buildMeasOp(index, '0');

auto tmp0 = conjugateTranspose(operation);
auto tmp1 = multiply(e, densityFromMatrixEdge(tmp0), 0, false);
auto tmp2 = multiply(densityFromMatrixEdge(operation), tmp1, 0, true);
auto densityMatrixTrace = trace(tmp2);

std::uniform_real_distribution<fp> dist(0., 1.);
if (const auto threshold = dist(mt); threshold > densityMatrixTrace.r) {
operation = buildMeasOp(index, '1');
tmp0 = conjugateTranspose(operation);
tmp1 = multiply(e, densityFromMatrixEdge(tmp0), 0, false);
tmp2 = multiply(densityFromMatrixEdge(operation), tmp1, 0, true);
measuredResult = '1';
densityMatrixTrace = trace(tmp2);
}

incRef(tmp2);
dEdge::alignDensityEdge(e);
decRef(e);
e = tmp2;
dEdge::setDensityMatrixTrue(e);

// Normalize density matrix
auto result = e.w / densityMatrixTrace;
cn.decRef(e.w);
e.w = cn.lookup(result);
cn.incRef(e.w);
return {e, measuredResult};
}

/**
* @brief Performs a specific measurement on the given state vector decision
* diagram. Collapses the state according to the measurement result.
Expand Down Expand Up @@ -1350,21 +1415,15 @@ template <class Config> class Package {
}
}

dEdge applyOperationToDensity(dEdge& e, const mEdge& operation,
const bool generateDensityMatrix = false) {
dEdge applyOperationToDensity(dEdge& e, const mEdge& operation) {
auto tmp0 = conjugateTranspose(operation);
auto tmp1 = multiply(e, densityFromMatrixEdge(tmp0), 0, false);
auto tmp2 = multiply(densityFromMatrixEdge(operation), tmp1, 0,
generateDensityMatrix);
auto tmp2 = multiply(densityFromMatrixEdge(operation), tmp1, 0, true);
incRef(tmp2);
dEdge::alignDensityEdge(e);
decRef(e);
e = tmp2;

if (generateDensityMatrix) {
dEdge::setDensityMatrixTrue(e);
}

dEdge::setDensityMatrixTrue(e);
return e;
}

Expand Down Expand Up @@ -1786,7 +1845,7 @@ template <class Config> class Package {
auto r = trace(a, eliminate);
return cn.lookup(r);
}
ComplexValue trace(const mEdge& a) {
template <class Edge> ComplexValue trace(const Edge& a) {
const auto eliminate = std::vector<bool>(nqubits, true);
return trace(a, eliminate).w;
}
Expand Down Expand Up @@ -1814,17 +1873,19 @@ template <class Config> class Package {

private:
/// TODO: introduce a compute table for the trace?
mCachedEdge trace(const mEdge& a, const std::vector<bool>& eliminate,
std::size_t alreadyEliminated = 0) {
template <class OperandNode>
CachedEdge<OperandNode> trace(const Edge<OperandNode>& a,
const std::vector<bool>& eliminate,
std::size_t alreadyEliminated = 0) {
const auto aWeight = static_cast<ComplexValue>(a.w);
if (aWeight.approximatelyZero()) {
return mCachedEdge::zero();
return CachedEdge<OperandNode>::zero();
}

if (mNode::isTerminal(a.p) ||
if (OperandNode::isTerminal(a.p) ||
std::none_of(eliminate.begin(), eliminate.end(),
[](bool v) { return v; })) {
return {a.p, aWeight};
return CachedEdge<OperandNode>{a.p, aWeight};
}

const auto v = a.p->v;
Expand All @@ -1837,12 +1898,12 @@ template <class Config> class Package {
return r;
}

std::array<mCachedEdge, NEDGE> edge{};
std::transform(
a.p->e.cbegin(), a.p->e.cend(), edge.begin(),
[this, &eliminate, &alreadyEliminated](const mEdge& e) -> mCachedEdge {
return trace(e, eliminate, alreadyEliminated);
});
std::array<CachedEdge<OperandNode>, NEDGE> edge{};
std::transform(a.p->e.cbegin(), a.p->e.cend(), edge.begin(),
[this, &eliminate, &alreadyEliminated](
const Edge<OperandNode>& e) -> CachedEdge<OperandNode> {
return trace(e, eliminate, alreadyEliminated);
});
const auto adjustedV =
static_cast<Qubit>(static_cast<std::size_t>(a.p->v) -
(static_cast<std::size_t>(std::count(
Expand Down
Loading
Loading