diff --git a/CMakeLists.txt b/CMakeLists.txt index 3349581b..56eaa95e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,14 +2,13 @@ cmake_minimum_required(VERSION 3.14...3.19) project(ddsim LANGUAGES CXX - VERSION 1.3.0 + VERSION 1.4.0 DESCRIPTION "DDSIM - A JKQ quantum simulator based on decision diagrams" ) set_property(GLOBAL PROPERTY USE_FOLDERS ON) option(COVERAGE "Configure for coverage report generation") option(GENERATE_POSITION_INDEPENDENT_CODE "Generate position independent code") -option(TEST_WITH_SANITIZERS "Run test with ASAN and UBSAN") set(default_build_type "Release") if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) diff --git a/LICENSE b/LICENSE index 1bb1f8bf..975d8e1b 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2020 Stefan Hillmich, Lukas Burgholzer, Robert Wille +Copyright (c) 2021 Stefan Hillmich, Lukas Burgholzer, Thomas Grurl, and Robert Wille Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal diff --git a/README.md b/README.md index bde58650..a6938304 100644 --- a/README.md +++ b/README.md @@ -10,8 +10,6 @@ A tool for quantum circuit simulation by the [Institute for Integrated Circuits](https://iic.jku.at/eda/) at the [Johannes Kepler University Linz](https://jku.at) and a part of the [JKQ toolset](https://github.com/iic-jku/jkq). -Developers: Stefan Hillmich, Lukas Burgholzer, Thomas Grurl, and Robert Wille. - The tool builds upon [our quantum functionality representation (QFR)](https://github.com/iic-jku/qfr.git) which in turns builds on [our decision diagram (DD) package](https://github.com/iic-jku/dd_package.git). For more information, on our work on quantum circuit simulation please visit [iic.jku.at/eda/research/quantum_simulation](https://iic.jku.at/eda/research/quantum_simulation) or, for more information on our work on noise-aware quantum circuit simulation, please visit [iic.jku.at/eda/research/noise_aware_simulation](https://iic.jku.at/eda/research/noise_aware_simulation). @@ -59,7 +57,7 @@ The following additional algorithms are integrated in [QFR](https://github.com/i For details on the available methods we refer to [iic.jku.at/eda/research/quantum_simulation](https://iic.jku.at/eda/research/quantum_simulation). -The simulator is based on [[1]](https://iic.jku.at/files/eda/2018_tcad_advanced_simulation_quantum_computation.pdf) and can either be used as a **standalone executable** with command-line interface, or as a **library** for the incorporation in other projects. +The simulator is based on the references listed below and can either be used as a **standalone executable** with command-line interface, or as a **library** for the incorporation in other projects. ## System Requirements @@ -114,29 +112,26 @@ The standalone executable is launched in the following way, showing available op ```commandline $ ./ddsim_simple --help JKQ DDSIM by https://iic.jku.at/eda/ -- Allowed options: --h [ --help ] produce help message ---seed arg (=0) seed for random number generator (default zero is possibly directly used as seed!) ---shots arg (=0) number of measurements (if the algorithm does not contain non-unitary gates, weak simulation is used) ---display_vector display the state vector ---ps print simulation stats (applied gates, sim. time, and maximal size of the DD) ---verbose Causes some simulators to print additional information to STDERR ---benchmark print simulation stats in a single CSV style line (overrides --ps and suppresses most other output, please don't rely on the format across versions) ---simulate_file arg simulate a quantum circuit given by file (detection by the file extension) ---simulate_qft arg simulate Quantum Fourier Transform for given number of qubits ---simulate_ghz arg simulate state preparation of GHZ state for given number of qubits ---step_fidelity arg (=1) target fidelity for each approximation run (>=1 = disable approximation) ---steps arg (=1) number of approximation steps ---initial_reorder arg (=0) Try to find a good initial variable order (0=None, 1=Most affected qubits to the top, 2=Most affected targets to the top) ---dynamic_reorder arg (=0) Apply reordering strategy during simulation (0=None, 1=Sifting, 2=Move2Top) ---post_reorder arg (=0) Apply a reordering strategy after simulation (0=None, 1=Sifting) ---simulate_grover arg simulate Grover's search for given number of qubits with random oracle ---simulate_grover_emulated arg simulate Grover's search for given number of qubits with random oracle and emulation ---simulate_grover_oracle_emulated arg simulate Grover's search for given number of qubits with given oracle and emulation ---simulate_shor arg simulate Shor's algorithm factoring this number ---simulate_shor_coprime arg (=0) coprime number to use with Shor's algorithm (zero randomly generates a coprime) ---simulate_shor_no_emulation Force Shor simulator to do modular exponentiation instead of using emulation (you'll usually want emulation) ---simulate_fast_shor arg simulate Shor's algorithm factoring this number with intermediate measurements ---simulate_fast_shor_coprime arg (=0) coprime number to use with Shor's algorithm (zero randomly generates a coprime) + -h [ --help ] produce help message + --seed arg (=0) seed for random number generator (default zero is possibly directly used as seed!) + --shots arg (=0) number of measurements (if the algorithm does not contain non-unitary gates, weak simulation is used) + --pv display the state vector as list of pairs (real and imaginary parts) + --ps print simulation stats (applied gates, sim. time, and maximal size of the DD) + --pm print measurement results + --verbose Causes some simulators to print additional information to STDERR + --simulate_file arg simulate a quantum circuit given by file (detection by the file extension) + --simulate_qft arg simulate Quantum Fourier Transform for given number of qubits + --simulate_ghz arg simulate state preparation of GHZ state for given number of qubits + --step_fidelity arg (=1) target fidelity for each approximation run (>=1 = disable approximation) + --steps arg (=1) number of approximation steps + --simulate_grover arg simulate Grover's search for given number of qubits with random oracle + --simulate_grover_emulated arg simulate Grover's search for given number of qubits with random oracle and emulation + --simulate_grover_oracle_emulated arg simulate Grover's search for given number of qubits with given oracle and emulation + --simulate_shor arg simulate Shor's algorithm factoring this number + --simulate_shor_coprime arg (=0) coprime number to use with Shor's algorithm (zero randomly generates a coprime) + --simulate_shor_no_emulation Force Shor simulator to do modular exponentiation instead of using emulation (you'll usually want emulation) + --simulate_fast_shor arg simulate Shor's algorithm factoring this number with intermediate measurements + --simulate_fast_shor_coprime arg (=0) coprime number to use with Shor's algorithm (zero randomly generates a coprime) ``` @@ -145,42 +140,25 @@ The output is JSON-formatted as shown below (with hopefully intuitive naming). ```commandline $ cmake -DCMAKE_BUILD_TYPE=Release -S . -B build $ cmake --build build --config Release --target ddsim_simple -$ ./build/ddsim_simple --simulate_ghz 4 --display_vector --shots 1000 --ps +$ ./build/ddsim_simple --simulate_ghz 4 --shots 1000 --ps --pm { - "measurements": { + "measurement_results": { "0000": 484, "1111": 516 }, - "state_vector": [ - +0.707107+0i, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - +0.707107+0i - ], - "non_zero_entries": 2, "statistics": { - "simulation_time": 0.000104, - "measurement_time": 0.000104, - "benchmark": "", - "shots": 1000, - "distinct_results": 2, - "n_qubits": 4, "applied_gates": 4, + "approximation_runs": "0", + "benchmark": "entanglement_4", + "distinct_results": 2, + "final_fidelity": "1.000000", "max_nodes": 9, - "path_of_least_resistance": "1111", - "seed": 0 + "n_qubits": 4, + "seed": "0", + "shots": 1000, + "simulation_time": 0.00013726699398830533, + "single_shots": "1", + "step_fidelity": "1.000000" } } ``` @@ -226,72 +204,50 @@ $ cmake --build build --config Release --target ddsim_noise_aware The simulator provides a help function which is called in the following way: ```commandline -$ ./build/ddsim_noise_aware --help +$ ./build/ddsim_noise_aware -h JKQ DDSIM by https://iic.jku.at/eda/ -- Allowed options: -h [ --help ] produce help message - --seed arg (=0) seed for random number generator - (default zero is possibly directly used - as seed!) - --shots arg (=0) number of measurements (if the - algorithm does not contain non-unitary - gates, weak simulation is used) - --display_vector display the state vector - --ps print simulation stats (applied gates, - sim. time, and maximal size of the DD) - --verbose Causes some simulators to print - additional information to STDERR - --benchmark print simulation stats in a single CSV - style line (overrides --ps and - suppresses most other output, please - don't rely on the format across - versions) - --simulate_file arg simulate a quantum circuit given by - file (detection by the file extension) - . - . - . - --noise_effects arg Noise effects (A (=amplitude damping),D - (=depolarization),P (=phase flip)) in - the form of a character string - describing the noise effects (default=" - ") - --noise_prob arg Probability for applying noise - (default=0.001) - --confidence arg Confidence in the error bound of the - stochastic simulation (default= 0.05) - --error_bound arg Error bound of the stochastic - simulation (default=0.1) - --stoch_runs arg (=0) Number of stochastic runs. When the - value is 0 the value is calculated - using the confidence, error_bound and - number of tracked properties. (default - = 0) - --properties arg Comma separated list of tracked - properties. Note that -1 is the - fidelity and "-" can be used to specify - a range. (default="0-1000") + --seed arg (=0) seed for random number generator (default zero is possibly directly used as seed!) + --pm print measurements + --ps print simulation stats (applied gates, sim. time, and maximal size of the DD) + --verbose Causes some simulators to print additional information to STDERR + --simulate_file arg simulate a quantum circuit given by file (detection by the file extension) + --step_fidelity arg (=1) target fidelity for each approximation run (>=1 = disable approximation) + --steps arg (=1) number of approximation steps + --noise_effects arg (=APD) Noise effects (A (=amplitude damping),D (=depolarization),P (=phase flip)) in the form of a character string describing the noise effects + --noise_prob arg (=0.001) Probability for applying noise + --confidence arg (=0.05) Confidence in the error bound of the stochastic simulation + --error_bound arg (=0.1) Error bound of the stochastic simulation + --stoch_runs arg (=0) Number of stochastic runs. When the value is 0 the value is calculated using the confidence, error_bound and number of tracked properties. + --properties arg (=-3-1000) Comma separated list of tracked properties. Note that -1 is the fidelity and "-" can be used to specify a range. + +Process finished with exit code 0 + ``` An example run, with amplitude damping, phase flip, and depolarization error (each with a probability of 0.1% whenever a gate is applied) looks like this: ```commandline -$ ./build/ddsim_noise_aware --noise_effects APD --stoch_runs 10000 --noise_prob 0.001 --simulate_file /home/user/adder4.qasm -Conducting perfect run ... -Conducting 10000 runs using 4 cores ... -Starting 4 threads -Calculating amplitudes from all runs ... -Probabilities are ... (probabilities < 0.001 are omitted) -state=|0000> proba=0.00975488 -state=|0001> proba=0.00784512 -state=|0100> proba=0.00175135 -state=|0101> proba=0.00194865 -state=|0110> proba=0.00404999 -state=|1000> proba=0.0229565 -state=|1001> proba=0.937243 -state=|1011> proba=0.00229939 -state=|1100> proba=0.00215185 -state=|1101> proba=0.00284815 -state=|1110> proba=0.00505003 +$ ./build/ddsim_noise_aware --ps --noise_effects APD --stoch_runs 10000 --noise_prob 0.001 --simulate_file adder4.qasm +{ + "statistics": { + "applied_gates": 23, + "approximation_runs": "0.000000", + "benchmark": "stoch_APD_adder_n4", + "final_fidelity": "0.937343", + "max_nodes": 10, + "mean_stoch_run_time": "0.015796", + "n_qubits": 4, + "parallel_instances": "28", + "perfect_run_time": "0.000066", + "seed": "0", + "simulation_time": 5.911194324493408, + "step_fidelity": "1.000000", + "stoch_runs": 10000, + "stoch_wall_time": "5.911118", + "threads": 28 + } +} ``` ## Running Tests @@ -329,7 +285,9 @@ If you use our tool for your research, we will be thankful if you refer to it by
-[1] A. Zulehner and R. Wille, “Advanced Simulation of Quantum Computations,” Transactions on CAD of Integrated Circuits and Systems (TCAD), vol. 38, no. 5, pp. 848–859, 2019 + + [1] A. Zulehner and R. Wille, “Advanced Simulation of Quantum Computations,” Transactions on CAD of Integrated Circuits and Systems (TCAD), vol. 38, no. 5, pp. 848–859, 2019 + ```bibtex @article{zulehner2019advanced, @@ -346,7 +304,29 @@ If you use our tool for your research, we will be thankful if you refer to it by
-[2] T. Grurl, R. Kueng, J. Fuß, and R. Wille, “Stochastic Quantum Circuit Simulation Using Decision Diagrams,” in Design, Automation and Test in Europe (DATE), 2021 + + [2] S. Hillmich, I. L. Markov, and R. Wille, “Just Like the Real Thing: Fast Weak Simulation of Quantum Computation,” in Design Automation Conference (DAC), 2020 + + +```bibtex +@inproceedings{DBLP:conf/dac/HillmichMW20, + author = {Stefan Hillmich and + Igor L. Markov and + Robert Wille}, + title = {Just Like the Real Thing: {F}ast Weak Simulation of Quantum Computation}, + booktitle = {Design Automation Conference}, + publisher = {{IEEE}}, + year = {2020} +} +``` +
+ + + +
+ + [3] T. Grurl, R. Kueng, J. Fuß, and R. Wille, “Stochastic Quantum Circuit Simulation Using Decision Diagrams,” in Design, Automation and Test in Europe (DATE), 2021 + ```bibtex @inproceedings{Grurl2020, @@ -361,4 +341,24 @@ If you use our tool for your research, we will be thankful if you refer to it by
+
+ + [4] S. Hillmich, R. Kueng, I. L. Markov, and R. Wille, "As Accurate as Needed, as Efficient as Possible: Approximations in DD-based Quantum Circuit Simulation," in Design, Automation and Test in Europe (DATE), 2021 + + +```bibtex +@inproceedings{DBLP:conf/date/HillmichKMW21, + author = {Stefan Hillmich and + Richard Kueng and + Igor L. Markov and + Robert Wille}, + title = {As Accurate as Needed, as Efficient as Possible: Approximations in DD-based Quantum Circuit Simulation}, + booktitle = {Design, Automation and Test in Europe}, + year = {2021} +} +``` +
+ + + diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index be14dc24..4b0cb805 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -20,10 +20,17 @@ macro(add_sim_executable APPNAME) set_target_properties(${PROJECT_NAME}_${APPNAME} PROPERTIES EXPORT_NAME ${PROJECT_NAME}_${APPNAME}) endmacro() -if (Boost_FOUND) +set(THREADS_PREFER_PTHREAD_FLAG ON) +find_package( Threads ) +link_libraries( Threads::Threads) + +if(Boost_FOUND) add_sim_executable(simple Boost::program_options) add_sim_executable(primebases Boost::program_options) - + if(Threads_FOUND) + add_sim_executable(noise_aware Boost::program_options ) + # target_link_libraries(${PROJECT_NAME}_noise_aware PUBLIC Threads::Threads) + endif() else () message(WARNING "Did not find Boost! Commandline interface will not be an available target!") endif () diff --git a/apps/noise_aware.cpp b/apps/noise_aware.cpp new file mode 100644 index 00000000..d3e8dff1 --- /dev/null +++ b/apps/noise_aware.cpp @@ -0,0 +1,110 @@ +#include "StochasticNoiseSimulator.hpp" +#include "nlohmann/json.hpp" + +#include +#include +#include +#include +#include +#include + +namespace nl = nlohmann; + +int main(int argc, char** argv) { + namespace po = boost::program_options; + unsigned long long seed; + + po::options_description description("JKQ DDSIM by https://iic.jku.at/eda/ -- Allowed options"); + // clang-format off + description.add_options() + ("help,h", "produce help message") + ("seed", po::value(&seed)->default_value(0), "seed for random number generator (default zero is possibly directly used as seed!)") + ("pm", "print measurements") + ("ps", "print simulation stats (applied gates, sim. time, and maximal size of the DD)") + ("verbose", "Causes some simulators to print additional information to STDERR") + + ("simulate_file", po::value(), "simulate a quantum circuit given by file (detection by the file extension)") + ("step_fidelity", po::value()->default_value(1.0), "target fidelity for each approximation run (>=1 = disable approximation)") + ("steps", po::value()->default_value(1), "number of approximation steps") + + ("noise_effects", po::value()->default_value("APD"), "Noise effects (A (=amplitude damping),D (=depolarization),P (=phase flip)) in the form of a character string describing the noise effects (default=\"APD\")") + ("noise_prob", po::value()->default_value(0.001), "Probability for applying noise (default=0.001)") + ("confidence", po::value()->default_value(0.05), "Confidence in the error bound of the stochastic simulation (default= 0.05)") + ("error_bound", po::value()->default_value(0.1), "Error bound of the stochastic simulation (default=0.1)") + ("stoch_runs", po::value()->default_value(0), "Number of stochastic runs. When the value is 0 the value is calculated using the confidence, error_bound and number of tracked properties. (default = 0)") + ("properties", po::value()->default_value("-3-1000"), R"(Comma separated list of tracked properties. Note that -1 is the fidelity and "-" can be used to specify a range. (default="-3-1000"))") + ; + // clang-format on + po::variables_map vm; + try { + po::store(po::parse_command_line(argc, argv, description), vm); + + if (vm.count("help")) { + std::cout << description; + return 0; + } + po::notify(vm); + } catch (const po::error& e) { + std::cerr << "[ERROR] " << e.what() << "! Try option '--help' for available commandline options.\n"; + std::exit(1); + } + + std::unique_ptr quantumComputation; + std::unique_ptr ddsim{nullptr}; + + if (vm.count("simulate_file")) { + const std::string fname = vm["simulate_file"].as(); + quantumComputation = std::make_unique(fname); + ddsim = std::make_unique(quantumComputation, + vm["steps"].as(), + vm["step_fidelity"].as(), + seed); + } else { + std::cerr << "Did not find anything to simulate. See help below.\n" + << description << "\n"; + std::exit(1); + } + + if (quantumComputation && quantumComputation->getNqubits() > 100) { + std::clog << "[WARNING] Quantum computation contains quite many qubits. You're jumping into the deep end.\n"; + } + + ddsim->setNoiseEffects(vm["noise_effects"].as()); + ddsim->setAmplitudeDampingProbability(vm["noise_prob"].as()); + ddsim->stoch_confidence = vm["confidence"].as(); + ddsim->setRecordedProperties(vm["properties"].as()); + ddsim->stoch_error_margin = vm["error_bound"].as(); + ddsim->stochastic_runs = vm["stoch_runs"].as(); + + auto t1 = std::chrono::steady_clock::now(); + + const std::map measurement_results = ddsim->StochSimulate(); + + auto t2 = std::chrono::steady_clock::now(); + + std::chrono::duration duration_simulation = t2 - t1; + + nl::json output_obj; + + if (vm.count("ps")) { + output_obj["statistics"] = { + {"simulation_time", duration_simulation.count()}, + {"benchmark", ddsim->getName()}, + {"stoch_runs", ddsim->stochastic_runs}, + {"threads", ddsim->max_instances}, + {"n_qubits", +ddsim->getNumberOfQubits()}, + {"applied_gates", ddsim->getNumberOfOps()}, + {"max_nodes", ddsim->getMaxNodeCount()}, + {"seed", ddsim->getSeed()}, + }; + + for (const auto& item: ddsim->AdditionalStatistics()) { + output_obj["statistics"][item.first] = item.second; + } + } + + if (vm.count("pm")) { + output_obj["measurement_results"] = measurement_results; + } + std::cout << std::setw(2) << output_obj << std::endl; +} diff --git a/apps/simple.cpp b/apps/simple.cpp index 4e0e9b1c..e00b556d 100644 --- a/apps/simple.cpp +++ b/apps/simple.cpp @@ -3,16 +3,19 @@ #include "ShorFastSimulator.hpp" #include "ShorSimulator.hpp" #include "Simulator.hpp" +#include "algorithms/Entanglement.hpp" +#include "algorithms/Grover.hpp" +#include "algorithms/QFT.hpp" +#include "nlohmann/json.hpp" -#include -#include -#include #include #include #include #include #include +namespace nl = nlohmann; + int main(int argc, char** argv) { namespace po = boost::program_options; // variables initialized by boost program_options default values @@ -24,30 +27,31 @@ int main(int argc, char** argv) { ApproximationInfo::ApproximationWhen approx_when; po::options_description description("JKQ DDSIM by https://iic.jku.at/eda/ -- Allowed options"); - description.add_options()("help,h", "produce help message")("seed", po::value<>(&seed)->default_value(0), - "seed for random number generator (default zero is possibly directly used as seed!)")("shots", po::value<>(&shots)->default_value(0), - "number of measurements (if the algorithm does not contain non-unitary gates, weak simulation is used)")("display_vector", "display the state vector")("ps", "print simulation stats (applied gates, sim. time, and maximal size of the DD)")("verbose", "Causes some simulators to print additional information to STDERR")("benchmark", - "print simulation stats in a single CSV style line (overrides --ps and suppresses most other output, please don't rely on the format across versions)") - - ("simulate_file", po::value(), - "simulate a quantum circuit given by file (detection by the file extension)")("simulate_qft", po::value(), "simulate Quantum Fourier Transform for given number of qubits")("simulate_ghz", po::value(), - "simulate state preparation of GHZ state for given number of qubits")("step_fidelity", po::value<>(&step_fidelity)->default_value(1.0), - "target fidelity for each approximation run (>=1 = disable approximation)")("steps", po::value<>(&approx_steps)->default_value(1), "number of approximation steps")("approx_when", po::value<>(&approx_when)->default_value(ApproximationInfo::FidelityDriven), - "approximation method ('fidelity' (default) or 'memory'")("approx_state", - "do excessive approximation runs at the end of the simulation to see how the quantum state behaves") - - ("simulate_grover", po::value(), - "simulate Grover's search for given number of qubits with random oracle")("simulate_grover_emulated", po::value(), - "simulate Grover's search for given number of qubits with random oracle and emulation")("simulate_grover_oracle_emulated", po::value(), - "simulate Grover's search for given number of qubits with given oracle and emulation") - - ("simulate_shor", po::value(), "simulate Shor's algorithm factoring this number")("simulate_shor_coprime", po::value()->default_value(0), - "coprime number to use with Shor's algorithm (zero randomly generates a coprime)")("simulate_shor_no_emulation", - "Force Shor simulator to do modular exponentiation instead of using emulation (you'll usually want emulation)") - - ("simulate_fast_shor", po::value(), - "simulate Shor's algorithm factoring this number with intermediate measurements")("simulate_fast_shor_coprime", po::value()->default_value(0), - "coprime number to use with Shor's algorithm (zero randomly generates a coprime)"); + // clang-format off + description.add_options() + ("help,h", "produce help message") + ("seed", po::value<>(&seed)->default_value(0), "seed for random number generator (default zero is possibly directly used as seed!)") + ("shots", po::value<>(&shots)->default_value(0), "number of measurements (if the algorithm does not contain non-unitary gates, weak simulation is used)") + ("pv", "display the state vector") + ("ps", "print simulation stats (applied gates, sim. time, and maximal size of the DD)") + ("pm", "print measurement results") + ("verbose", "Causes some simulators to print additional information to STDERR") + ("simulate_file", po::value(), "simulate a quantum circuit given by file (detection by the file extension)") + ("simulate_qft", po::value(), "simulate Quantum Fourier Transform for given number of qubits") + ("simulate_ghz", po::value(), "simulate state preparation of GHZ state for given number of qubits") + ("step_fidelity", po::value<>(&step_fidelity)->default_value(1.0), "target fidelity for each approximation run (>=1 = disable approximation)") + ("steps", po::value<>(&approx_steps)->default_value(1), "number of approximation steps") + ("approx_when", po::value<>(&approx_when)->default_value(ApproximationInfo::FidelityDriven), "approximation method ('fidelity' (default) or 'memory'") + ("approx_state", "do excessive approximation runs at the end of the simulation to see how the quantum state behaves") + ("simulate_grover", po::value(), "simulate Grover's search for given number of qubits with random oracle") + ("simulate_grover_emulated", po::value(), "simulate Grover's search for given number of qubits with random oracle and emulation") + ("simulate_grover_oracle_emulated", po::value(), "simulate Grover's search for given number of qubits with given oracle and emulation") + ("simulate_shor", po::value(), "simulate Shor's algorithm factoring this number") + ("simulate_shor_coprime", po::value()->default_value(0), "coprime number to use with Shor's algorithm (zero randomly generates a coprime)") + ("simulate_shor_no_emulation", "Force Shor simulator to do modular exponentiation instead of using emulation (you'll usually want emulation)") + ("simulate_fast_shor", po::value(), "simulate Shor's algorithm factoring this number with intermediate measurements") + ("simulate_fast_shor_coprime", po::value()->default_value(0),"coprime number to use with Shor's algorithm (zero randomly generates a coprime)"); + // clang-format on po::variables_map vm; try { po::store(po::parse_command_line(argc, argv, description), vm); @@ -184,69 +188,32 @@ int main(int argc, char** argv) { } } - if (vm.count("benchmark")) { - auto more_info = ddsim->AdditionalStatistics(); - std::cout << ddsim->getName() << ", " - << ddsim->getNumberOfQubits() << ", " - << std::fixed << duration_simulation.count() << std::defaultfloat << ", " - << more_info["single_shots"] << "," - << more_info["approximation_runs"] << "," - << more_info["final_fidelity"] << ", " - << more_info["coprime_a"] << ", " - << more_info["sim_result"] << ", " - << more_info["polr_result"] << ", " - << ddsim->getSeed() << ", " - << ddsim->getNumberOfOps() << ", " - << ddsim->getMaxNodeCount() - << "\n"; - return 0; - } - - std::cout << "{\n"; + nl::json output_obj; - if (!m.empty()) { - std::cout << " \"measurements\": {"; - bool first_element = true; - for (const auto& element: m) { - std::cout << (first_element ? "" : ",") << "\n \"" << element.first << "\": " << element.second; - first_element = false; - } - std::cout << "\n },\n"; + if (vm.count("pm")) { + output_obj["measurement_results"] = m; } - if (vm.count("display_vector")) { - std::cout << " \"state_vector\": ["; - - bool first_element = true; - unsigned long long non_zero_entries = 0; - for (const auto& element: ddsim->getVector()) { - if (element.r != 0 || element.i != 0) { - non_zero_entries++; - std::cout << (first_element ? "" : ",") << "\n " << std::showpos << element.r << element.i << "i" - << std::noshowpos; - } else { - std::cout << (first_element ? "" : ",") << "\n 0"; - } - first_element = false; - } - std::cout << "\n ],\n"; - std::cout << " \"non_zero_entries\": " << non_zero_entries << ",\n"; + + if (vm.count("pv")) { + output_obj["state_vector"] = ddsim->getVectorPair(); } if (vm.count("ps")) { - std::cout << " \"statistics\": {\n" - << " \"simulation_time\": " << std::fixed << duration_simulation.count() << std::defaultfloat - << ",\n" - << " \"benchmark\": \"" << ddsim->getName() << "\",\n" - << " \"shots\": " << shots << ",\n" - << " \"distinct_results\": " << m.size() << ",\n" - << " \"n_qubits\": " << ddsim->getNumberOfQubits() << ",\n" - << " \"applied_gates\": " << ddsim->getNumberOfOps() << ",\n" - << " \"max_nodes\": " << ddsim->getMaxNodeCount() << ",\n"; + output_obj["statistics"] = { + {"simulation_time", duration_simulation.count()}, + {"benchmark", ddsim->getName()}, + {"n_qubits", +ddsim->getNumberOfQubits()}, + {"applied_gates", ddsim->getNumberOfOps()}, + {"max_nodes", ddsim->getMaxNodeCount()}, + {"shots", shots}, + {"distinct_results", m.size()}, + {"seed", ddsim->getSeed()}, + }; + for (const auto& item: ddsim->AdditionalStatistics()) { - std::cout << " \"" << item.first << "\": \"" << item.second << "\",\n"; + output_obj["statistics"][item.first] = item.second; } - std::cout << " \"seed\": " << ddsim->getSeed() << "\n" - << " },\n"; } - std::cout << " \"dummy\": 0\n}\n"; // trailing element to make json printout easier + + std::cout << std::setw(2) << output_obj << std::endl; } diff --git a/extern/qfr b/extern/qfr index d420683f..7bb2681a 160000 --- a/extern/qfr +++ b/extern/qfr @@ -1 +1 @@ -Subproject commit d420683fb9e18a522d60d318f6e9311886acd051 +Subproject commit 7bb2681a3fcecd42f7440fb352ce86900cf0fcb4 diff --git a/include/Simulator.hpp b/include/Simulator.hpp index 327daeb3..74cba7bc 100644 --- a/include/Simulator.hpp +++ b/include/Simulator.hpp @@ -45,6 +45,8 @@ class Simulator { [[nodiscard]] std::vector getVector() const; + [[nodiscard]] std::vector> getVectorPair() const; + [[nodiscard]] std::size_t getActiveNodeCount() const { return dd->vUniqueTable.getActiveNodeCount(); } [[nodiscard]] std::size_t getMaxNodeCount() const { return dd->vUniqueTable.getMaxActiveNodes(); } diff --git a/include/StochasticNoiseSimulator.hpp b/include/StochasticNoiseSimulator.hpp new file mode 100644 index 00000000..87ff3324 --- /dev/null +++ b/include/StochasticNoiseSimulator.hpp @@ -0,0 +1,142 @@ +#ifndef DDSIM_STOCHASTICNOISESIMULATOR_HPP +#define DDSIM_STOCHASTICNOISESIMULATOR_HPP + +#include "QuantumComputation.hpp" +#include "Simulator.hpp" + +#include +#include +#include +#include +#include +#include +#include + +class StochasticNoiseSimulator: public Simulator { +public: + StochasticNoiseSimulator(std::unique_ptr& qc, const unsigned int step_number, const double step_fidelity): + qc(qc), step_number(step_number), step_fidelity(step_fidelity) { + dd->resize(qc->getNqubits()); + } + + StochasticNoiseSimulator(std::unique_ptr& qc, const unsigned int step_number, const double step_fidelity, unsigned long long seed): + Simulator(seed), qc(qc), step_number(step_number), step_fidelity(step_fidelity) { + dd->resize(qc->getNqubits()); + } + + StochasticNoiseSimulator(std::unique_ptr& qc, + const std::string& noise_effects, + double noise_prob, + long stoch_runs, + unsigned int step_number, + double step_fidelity, + const std::string& recorded_properties): + qc(qc), + step_number(step_number), step_fidelity(step_fidelity) { + setNoiseEffects(noise_effects); + setRecordedProperties(recorded_properties); + setAmplitudeDampingProbability(noise_prob); + stochastic_runs = stoch_runs; + } + + std::map Simulate(unsigned int shots) override; + + std::map StochSimulate(); + + std::map AdditionalStatistics() override { + return { + {"step_fidelity", std::to_string(step_fidelity)}, + {"approximation_runs", std::to_string(approximation_runs)}, + {"final_fidelity", std::to_string(final_fidelity)}, + {"perfect_run_time", std::to_string(perfect_run_time)}, + {"stoch_wall_time", std::to_string(stoch_run_time)}, + {"mean_stoch_run_time", std::to_string(mean_stoch_time)}, + {"parallel_instances", std::to_string(max_instances)}, + }; + }; + + [[nodiscard]] dd::QubitCount getNumberOfQubits() const override { return qc->getNqubits(); }; + + [[nodiscard]] std::size_t getNumberOfOps() const override { return qc->getNops(); }; + + [[nodiscard]] std::string getName() const override { return "stoch_" + gate_noise_types + "_" + qc->getName(); }; + + dd::Package::vEdge MeasureOneCollapsingConcurrent(unsigned short index, const std::unique_ptr& localDD, + dd::Package::vEdge local_root_edge, + std::mt19937_64& generator, + char* result, + bool assume_probability_normalization = true); + + void setNoiseEffects(const std::string& cGateNoise) { gate_noise_types = cGateNoise; } + + double noise_probability = 0.0; + dd::ComplexValue sqrt_amplitude_damping_probability{}; + dd::ComplexValue one_minus_sqrt_amplitude_damping_probability{}; + + void setAmplitudeDampingProbability(double cGateNoiseProbability) { + //The probability of amplitude damping (t1) often is double the probability , of phase flip, which is why I double it here + noise_probability = cGateNoiseProbability; + sqrt_amplitude_damping_probability = {sqrt(noise_probability * 2), 0}; + one_minus_sqrt_amplitude_damping_probability = {sqrt(1 - noise_probability * 2), 0}; + } + + double stoch_error_margin = 0.01; + double stoch_confidence = 0.05; + unsigned int stochastic_runs = 0; + + void setRecordedProperties(const std::string& input); + + std::vector> recorded_properties; + std::vector> recorded_properties_per_instance; + + std::string gate_noise_types; + + const unsigned int max_instances = std::max(1, static_cast(std::thread::hardware_concurrency()) - 4); + +private: + std::unique_ptr& qc; + + const unsigned int step_number; + const double step_fidelity; + double approximation_runs{0}; + long double final_fidelity{1.0L}; + + float perfect_run_time{0}; + float stoch_run_time{0}; + double mean_stoch_time{0}; + + void perfect_simulation_run(); + + void runStochSimulationForId(unsigned int stochRun, + int n_qubits, + dd::Package::vEdge rootEdgePerfectRun, + std::vector& recordedPropertiesStorage, + std::vector>& recordedPropertiesList, + unsigned long long local_seed); + + dd::Package::mEdge generateNoiseOperation(bool amplitudeDamping, + dd::Qubit target, + std::mt19937_64& engine, + std::uniform_real_distribution& distribution, + dd::Package::mEdge dd_operation, + std::unique_ptr& localDD); + + void applyNoiseOperation(const qc::Targets& targets, + const dd::Controls& control_qubits, + dd::Package::mEdge dd_op, + std::unique_ptr& localDD, + dd::Package::vEdge& local_root_edge, + std::mt19937_64& generator, + std::uniform_real_distribution& dist, + dd::Package::mEdge identityDD); + + [[nodiscard]] dd::NoiseOperationKind ReturnNoiseOperation(char i, double d) const; + + [[nodiscard]] std::string intToString(long target_number) const; + + double ApproximateEdgeByFidelity(std::unique_ptr& localDD, dd::Package::vEdge& edge, double targetFidelity, bool allLevels, bool removeNodes); + + dd::Package::vEdge RemoveNodesInPackage(std::unique_ptr& localDD, dd::Package::vEdge e, std::map& dag_edges); +}; + +#endif //DDSIM_STOCHASTICNOISESIMULATOR_HPP diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 54134b9b..9af78af0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -9,6 +9,8 @@ add_library(${PROJECT_NAME} ${CMAKE_CURRENT_SOURCE_DIR}/ShorSimulator.cpp ${PROJECT_SOURCE_DIR}/include/ShorFastSimulator.hpp ${CMAKE_CURRENT_SOURCE_DIR}/ShorFastSimulator.cpp + ${PROJECT_SOURCE_DIR}/include/StochasticNoiseSimulator.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/StochasticNoiseSimulator.cpp ) target_include_directories(${PROJECT_NAME} PUBLIC $) # set required C++ standard and disable compiler specific extensions diff --git a/src/Simulator.cpp b/src/Simulator.cpp index 321db728..cb67b6d5 100644 --- a/src/Simulator.cpp +++ b/src/Simulator.cpp @@ -1,5 +1,6 @@ #include "Simulator.hpp" +#include #include #include #include @@ -88,6 +89,20 @@ std::vector Simulator::getVector() const { return results; } +std::vector> Simulator::getVectorPair() const { + assert(getNumberOfQubits() < 60); // On 64bit system the vector can hold up to (2^60)-1 elements, if memory permits + std::string path(getNumberOfQubits(), '0'); + std::vector> results{1ull << getNumberOfQubits()}; + + for (unsigned long long i = 0; i < 1ull << getNumberOfQubits(); ++i) { + const std::string corrected_path{path.rbegin(), path.rend()}; + const dd::ComplexValue cv = dd->getValueByPath(root_edge, corrected_path); + results[i] = std::make_pair(cv.r, cv.i); + NextPath(path); + } + return results; +} + void Simulator::NextPath(std::string& s) { std::string::reverse_iterator iter = s.rbegin(), end = s.rend(); int carry = 1; @@ -209,12 +224,12 @@ char Simulator::MeasureOneCollapsing(const dd::Qubit index, const bool assume_pr dd::Package::mEdge m_gate = dd->makeGateDD(measure_m, getNumberOfQubits(), index); dd::Package::vEdge e = dd->multiply(m_gate, root_edge); - dd->decRef(root_edge); dd::Complex c = dd->cn.getTemporary(std::sqrt(1.0 / norm_factor), 0); CN::mul(c, e.w, c); e.w = dd->cn.lookup(c); dd->incRef(e); + dd->decRef(root_edge); root_edge = e; return result; diff --git a/src/StochasticNoiseSimulator.cpp b/src/StochasticNoiseSimulator.cpp new file mode 100644 index 00000000..6faa8f38 --- /dev/null +++ b/src/StochasticNoiseSimulator.cpp @@ -0,0 +1,727 @@ +#include "StochasticNoiseSimulator.hpp" + +#include +#include + +std::map StochasticNoiseSimulator::Simulate(unsigned int shots) { + bool has_nonunitary = false; + for (auto& op: *qc) { + if (op->isNonUnitaryOperation()) { + has_nonunitary = true; + break; + } + } + + if (!has_nonunitary) { + perfect_simulation_run(); + return MeasureAllNonCollapsing(shots); + } + + std::map m_counter; + + for (unsigned int i = 0; i < shots; i++) { + perfect_simulation_run(); + m_counter[MeasureAll()]++; + } + + return m_counter; +} + +void StochasticNoiseSimulator::perfect_simulation_run() { + const unsigned short n_qubits = qc->getNqubits(); + + root_edge = dd->makeZeroState(n_qubits); + dd->incRef(root_edge); + + unsigned long op_num = 0; + std::map classic_values; + + for (auto& op: *qc) { + if (op->isNonUnitaryOperation()) { + if (auto* nu_op = dynamic_cast(op.get())) { + if (op->getType() == qc::Measure) { + auto quantum = nu_op->getTargets(); + auto classic = nu_op->getClassics(); + + assert(quantum.size() == classic.size()); // this should not happen do to check in Simulate + + for (unsigned int i = 0; i < quantum.size(); ++i) { + auto result = MeasureOneCollapsing(quantum.at(i)); + assert(result == '0' || result == '1'); + classic_values[classic.at(i)] = (result == '1'); + } + } else if (op->getType() == qc::Barrier) { + continue; + } else { + throw std::runtime_error(std::string("Unsupported non-unitary functionality '") + op->getName() + "'."); + } + } else { + throw std::runtime_error(std::string("Dynamic cast to NonUnitaryOperation failed for '") + op->getName() + "'."); + } + dd->garbageCollect(); + } else { + if (op->isClassicControlledOperation()) { + if (auto* cc_op = dynamic_cast(op.get())) { + const auto start_index = static_cast(cc_op->getParameter().at(0)); + const auto length = static_cast(cc_op->getParameter().at(1)); + const unsigned int expected_value = cc_op->getExpectedValue(); + + unsigned int actual_value = 0; + for (unsigned int i = 0; i < length; i++) { + actual_value |= (classic_values[start_index + i] ? 1u : 0u) << i; + } + + //std::clog << "expected " << expected_value << " and actual value was " << actual_value << "\n"; + + if (actual_value != expected_value) { + continue; + } + } else { + throw std::runtime_error("Dynamic cast to ClassicControlledOperation failed."); + } + } + + auto dd_op = op->getDD(dd); + auto tmp = dd->multiply(dd_op, root_edge); + dd->incRef(tmp); + dd->decRef(root_edge); + root_edge = tmp; + + dd->garbageCollect(); + } + op_num++; + } +} + +std::map StochasticNoiseSimulator::StochSimulate() { + const unsigned short n_qubits = qc->getNqubits(); + + //Ceiling[(Log[estimatesProp] + Log[(2/confidence)])/(2*errorBound^2)] + if (stochastic_runs == 0) { + stochastic_runs = std::ceil((std::log(recorded_properties.size()) + std::log(2 / stoch_confidence)) / (2 * stoch_error_margin * stoch_error_margin)); + } + //std::clog << "Conducting perfect run...\n"; + const auto t1_perfect = std::chrono::steady_clock::now(); + perfect_simulation_run(); + const auto t2_perfect = std::chrono::steady_clock::now(); + perfect_run_time = std::chrono::duration(t2_perfect - t1_perfect).count(); + + dd::Edge noiseless_root_edge = root_edge; + + // Generate a vector for each instance + 1. the final vector stores the average of all runs and is calculated + // after the runs have finished + for (unsigned int i = 0; i < max_instances + 1; i++) { + recorded_properties_per_instance.emplace_back(std::vector(recorded_properties.size(), 0)); + } + //std::clog << "Conducting " << stochastic_runs << " runs using " << max_instances << " cores...\n"; + std::vector threadArray; + // The stochastic runs are applied in parallel + //std::clog << "Starting " << max_instances << " threads\n"; + const auto t1_stoch = std::chrono::steady_clock::now(); + for (unsigned int runID = 0; runID < max_instances; runID++) { + threadArray.emplace_back(&StochasticNoiseSimulator::runStochSimulationForId, + this, + runID, + n_qubits, + noiseless_root_edge, + std::ref(recorded_properties_per_instance[runID]), + std::ref(recorded_properties), + static_cast(mt())); + } + // wait for threads to finish + for (auto& thread: threadArray) { + thread.join(); + } + const auto t2_stoch = std::chrono::steady_clock::now(); + stoch_run_time = std::chrono::duration(t2_stoch - t1_stoch).count(); + //std::clog <<"Calculating amplitudes from all runs...\n"; + for (unsigned long j = 0; j < recorded_properties.size(); j++) { + for (unsigned int i = 0; i < max_instances; i++) { + recorded_properties_per_instance[max_instances][j] += recorded_properties_per_instance[i][j]; + } + recorded_properties_per_instance[max_instances][j] /= stochastic_runs; + } + + //std::clog << "Probabilities are ... (probabilities < 0.001 are omitted)\n"; + std::map measure_result; + for (unsigned long m = 0; m < recorded_properties.size(); m++) { + if (std::get<0>(recorded_properties[m]) == -3) { + mean_stoch_time = recorded_properties_per_instance[max_instances][m]; + } else if (std::get<0>(recorded_properties[m]) == -2) { + approximation_runs = recorded_properties_per_instance[max_instances][m]; + } else if (std::get<0>(recorded_properties[m]) == -1) { + final_fidelity = recorded_properties_per_instance[max_instances][m]; + } else if (recorded_properties_per_instance[max_instances][m] > 0.001 || m < 2) { + std::string amplitude = std::get<1>(recorded_properties[m]); + std::replace(amplitude.begin(), amplitude.end(), '2', '1'); + measure_result.insert({amplitude, recorded_properties_per_instance[max_instances][m]}); + } + } + return measure_result; +} + +dd::Package::vEdge StochasticNoiseSimulator::RemoveNodesInPackage(std::unique_ptr& localDD, dd::Package::vEdge e, std::map& dag_edges) { + if (e.isTerminal()) { + return e; + } + + const auto it = dag_edges.find(e.p); + if (it != dag_edges.end()) { + dd::Package::vEdge r = it->second; + if (r.w.approximatelyZero()) { + return dd::Package::vEdge::zero; + } + dd::Complex c = localDD->cn.getTemporary(); + dd::ComplexNumbers::mul(c, e.w, r.w); + r.w = localDD->cn.lookup(c); + return r; + } + + std::array edges{ + RemoveNodesInPackage(localDD, e.p->e.at(0), dag_edges), + RemoveNodesInPackage(localDD, e.p->e.at(1), dag_edges)}; + + dd::Package::vEdge r = localDD->makeDDNode(e.p->v, edges, false); + dag_edges[e.p] = r; + dd::Complex c = localDD->cn.getCached(); + dd::ComplexNumbers::mul(c, e.w, r.w); + r.w = localDD->cn.lookup(c); + return r; +} + +double StochasticNoiseSimulator::ApproximateEdgeByFidelity(std::unique_ptr& localDD, dd::Package::vEdge& edge, double targetFidelity, bool allLevels, bool removeNodes) { + std::queue q; + std::map probsMone; + + probsMone[edge.p] = dd::ComplexNumbers::mag2(edge.w); + + q.push(edge.p); + + while (!q.empty()) { + dd::Package::vNode* ptr = q.front(); + q.pop(); + const dd::fp parent_prob = probsMone[ptr]; + + if (ptr->e.at(0).w != dd::Complex::zero) { + if (probsMone.find(ptr->e.at(0).p) == probsMone.end()) { + q.push(ptr->e.at(0).p); + probsMone[ptr->e.at(0).p] = 0; + } + probsMone[ptr->e.at(0).p] = probsMone.at(ptr->e.at(0).p) + parent_prob * dd::ComplexNumbers::mag2(ptr->e.at(0).w); + } + + if (ptr->e.at(1).w != dd::Complex::zero) { + if (probsMone.find(ptr->e.at(1).p) == probsMone.end()) { + q.push(ptr->e.at(1).p); + probsMone[ptr->e.at(1).p] = 0; + } + probsMone[ptr->e.at(1).p] = probsMone.at(ptr->e.at(1).p) + parent_prob * dd::ComplexNumbers::mag2(ptr->e.at(1).w); + } + } + + std::vector nodes(getNumberOfQubits(), 0); + + std::vector, std::vector>>> qq( + getNumberOfQubits()); + + for (auto& it: probsMone) { + if (it.first->v < 0) { + continue; // ignore the terminal node which has v == -1 + } + nodes.at(it.first->v)++; + qq.at(it.first->v).emplace(1 - it.second, it.first); + } + + probsMone.clear(); + std::vector nodes_to_remove; + + int max_remove = 0; + for (int i = 0; i < getNumberOfQubits(); i++) { + double sum = 0.0; + int remove = 0; + std::vector tmp; + + while (!qq.at(i).empty()) { + auto p = qq.at(i).top(); + qq.at(i).pop(); + sum += 1 - p.first; + if (sum < 1 - targetFidelity) { + remove++; + if (allLevels) { + nodes_to_remove.push_back(p.second); + } else { + tmp.push_back(p.second); + } + } else { + break; + } + } + if (!allLevels && remove * i > max_remove) { + max_remove = remove * i; + nodes_to_remove = tmp; + } + } + + std::map dag_edges; + for (auto& it: nodes_to_remove) { + dag_edges[it] = dd::Package::vEdge::zero; + } + + dd::Package::vEdge newEdge = RemoveNodesInPackage(localDD, edge, dag_edges); + assert(!std::isnan(dd::CTEntry::val(edge.w.r))); + assert(!std::isnan(dd::CTEntry::val(edge.w.i))); + dd::Complex c = localDD->cn.getCached(std::sqrt(dd::ComplexNumbers::mag2(newEdge.w)), 0); + dd::ComplexNumbers::div(c, newEdge.w, c); + newEdge.w = localDD->cn.lookup(c); + + dd::fp fidelity = 0; + if (edge.p->v == newEdge.p->v) { + fidelity = localDD->fidelity(edge, newEdge); + } + + if (removeNodes) { + localDD->decRef(edge); + localDD->incRef(newEdge); + edge = newEdge; + } + return fidelity; +} + +void StochasticNoiseSimulator::runStochSimulationForId(unsigned int stochRun, + int n_qubits, + dd::Package::vEdge rootEdgePerfectRun, + std::vector& recordedPropertiesStorage, + std::vector>& recordedPropertiesList, + unsigned long long local_seed) { + std::mt19937_64 generator(local_seed); + std::uniform_real_distribution dist(0.0, 1.0L); + + const unsigned long numberOfRuns = stochastic_runs / max_instances + (stochRun < stochastic_runs % max_instances ? 1 : 0); + const int approx_mod = std::ceil(static_cast(qc->getNops()) / (step_number + 1)); + + //printf("Running %d times and using the dd at %p, using the cn object at %p\n", numberOfRuns, (void *) &localDD, (void *) &localDD->cn); + for (unsigned long current_run = 0; current_run < numberOfRuns; current_run++) { + const auto t1 = std::chrono::steady_clock::now(); + + std::unique_ptr localDD = std::make_unique(getNumberOfQubits()); + //dd::NoiseOperationTable noiseOperationTable(getNumberOfQubits()); + + dd::Package::mEdge identity_DD = localDD->makeIdent(n_qubits); + localDD->incRef(identity_DD); + dd::Package::vEdge local_root_edge = localDD->makeZeroState(n_qubits); + localDD->incRef(local_root_edge); + + std::map classic_values; + + unsigned int op_count = 0; + unsigned int approx_count = 0; + + for (auto& op: *qc) { + if (!op->isUnitary() && !(op->isClassicControlledOperation())) { + if (auto* nu_op = dynamic_cast(op.get())) { + if (nu_op->getName()[0] == 'M') { + auto quantum = nu_op->getTargets(); + auto classic = nu_op->getClassics(); + + assert(quantum.size() == classic.size()); // this should not happen do to check in Simulate + + for (unsigned int i = 0; i < quantum.size(); ++i) { + char result; + local_root_edge = MeasureOneCollapsingConcurrent(quantum.at(i), localDD, local_root_edge, generator, &result); + assert(result == '0' || result == '1'); + classic_values[classic.at(i)] = (result == '1'); + } + } else if (nu_op->getName()[0] == 'R' && nu_op->getName()[1] == 's' && nu_op->getName()[2] == 't') { + // Reset qubit + throw std::runtime_error("Warning: Reset is currently not supported"); + } else { + //Skipping barrier + if (op->getType() == qc::Barrier) { + continue; + } + throw std::runtime_error("Unsupported non-unitary functionality."); + } + } else { + throw std::runtime_error("Dynamic cast to NonUnitaryOperation failed."); + } + localDD->garbageCollect(false); + } else { + dd::Package::mEdge dd_op{}; + qc::Targets targets; + dd::Controls controls; + if (op->isClassicControlledOperation()) { + // Check if the operation is controlled by a classical register + auto* classic_op = dynamic_cast(op.get()); + bool execute_op = true; + auto expValue = classic_op->getExpectedValue(); + + for (unsigned int i = classic_op->getControlRegister().first; i < classic_op->getControlRegister().second; i++) { + if (classic_values[i] != (expValue % 2)) { + execute_op = false; + break; + } else { + expValue = expValue >> 1u; + } + } + dd_op = classic_op->getOperation()->getDD(localDD); + targets = classic_op->getOperation()->getTargets(); + controls = classic_op->getOperation()->getControls(); + if (!execute_op) { + // applyNoiseOperation(targets.front(), controls, identity_DD, (std::unique_ptr &) localDD, local_root_edge, generator, dist, line, identity_DD); + continue; + } + } else { + dd_op = op->getDD(localDD); + targets = op->getTargets(); + controls = op->getControls(); + } + + applyNoiseOperation(targets, controls, dd_op, localDD, local_root_edge, generator, dist, identity_DD); + if (step_fidelity < 1 && (op_count + 1) % approx_mod == 0) { + ApproximateEdgeByFidelity(localDD, local_root_edge, step_fidelity, false, true); + approx_count++; + } + //localDD->garbageCollect(false); + } + op_count++; + } + localDD->decRef(local_root_edge); + const auto t2 = std::chrono::steady_clock::now(); + + for (unsigned long i = 0; i < recordedPropertiesStorage.size(); i++) { + if (std::get<0>(recordedPropertiesList[i]) == -3) { + recordedPropertiesStorage[i] += std::chrono::duration(t2 - t1).count(); + } else if (std::get<0>(recordedPropertiesList[i]) == -2) { + recordedPropertiesStorage[i] += approx_count; + } else if (std::get<0>(recordedPropertiesList[i]) == -1) { + recordedPropertiesStorage[i] += localDD->fidelity(local_root_edge, rootEdgePerfectRun); + } else { + // extract amplitude for state + const auto basisVector = std::get<1>(recordedPropertiesList[i]); + auto amplitude = localDD->getValueByPath(local_root_edge, basisVector); + auto prob = amplitude.r * amplitude.r + amplitude.i * amplitude.i; + recordedPropertiesStorage[i] += prob; + } + } + } +} + +void StochasticNoiseSimulator::applyNoiseOperation(const qc::Targets& targets, + const dd::Controls& control_qubits, + dd::Package::mEdge dd_op, + std::unique_ptr& localDD, + dd::Package::vEdge& local_root_edge, + std::mt19937_64& generator, + std::uniform_real_distribution& dist, + dd::Package::mEdge identityDD) { + // TODO: Is this a meaningful way to handle multi-target operations? + for (auto& target: targets) { + auto operation = generateNoiseOperation(false, target, generator, dist, dd_op, localDD); + auto tmp = localDD->multiply(operation, local_root_edge); + + if (dd::ComplexNumbers::mag2(tmp.w) < dist(generator)) { + operation = generateNoiseOperation(true, target, generator, dist, dd_op, localDD); + tmp = localDD->multiply(operation, local_root_edge); + } + + if (tmp.w != dd::Complex::one) { + tmp.w = dd::Complex::one; + } + localDD->incRef(tmp); + localDD->decRef(local_root_edge); + local_root_edge = tmp; + } + // Apply noise to control qubits + for (auto control: control_qubits) { + auto operation = generateNoiseOperation(false, control.qubit, generator, dist, identityDD, localDD); + auto tmp = localDD->multiply(operation, local_root_edge); + if (dd::ComplexNumbers::mag2(tmp.w) < dist(generator)) { + operation = generateNoiseOperation(true, control.qubit, generator, dist, identityDD, localDD); + tmp = localDD->multiply(operation, local_root_edge); + } + if (tmp.w != dd::Complex::one) { + tmp.w = dd::Complex::one; + } + localDD->incRef(tmp); + localDD->decRef(local_root_edge); + local_root_edge = tmp; + } +} + +dd::Package::mEdge StochasticNoiseSimulator::generateNoiseOperation(bool amplitudeDamping, + dd::Qubit target, + std::mt19937_64& engine, + std::uniform_real_distribution& distribution, + dd::Package::mEdge dd_operation, + std::unique_ptr& localDD) { + dd::NoiseOperationKind effect; + + for (const auto& noise_type: gate_noise_types) { + if (noise_type != 'A') { + effect = ReturnNoiseOperation(noise_type, distribution(engine)); + } else { + if (amplitudeDamping) { + effect = dd::ATrue; + } else { + effect = dd::AFalse; + } + } + if (effect != dd::I) { + switch (effect) { + case (dd::ATrue): { + auto tmp_op = localDD->noiseOperationTable.lookup(getNumberOfQubits(), dd::ATrue, target); + if (tmp_op.p == nullptr) { + tmp_op = localDD->makeGateDD(dd::GateMatrix({dd::complex_zero, sqrt_amplitude_damping_probability, dd::complex_zero, dd::complex_zero}), getNumberOfQubits(), target); + localDD->noiseOperationTable.insert(dd::ATrue, target, tmp_op); + } + dd_operation = localDD->multiply(tmp_op, dd_operation); + break; + } + case (dd::AFalse): { + auto tmp_op = localDD->noiseOperationTable.lookup(getNumberOfQubits(), dd::AFalse, target); + if (tmp_op.p == nullptr) { + tmp_op = localDD->makeGateDD(dd::GateMatrix({dd::complex_one, dd::complex_zero, dd::complex_zero, one_minus_sqrt_amplitude_damping_probability}), getNumberOfQubits(), target); + localDD->noiseOperationTable.insert(dd::AFalse, target, tmp_op); + } + dd_operation = localDD->multiply(tmp_op, dd_operation); + break; + } + case (dd::X): { + dd_operation = localDD->multiply(localDD->makeGateDD(dd::Xmat, getNumberOfQubits(), target), dd_operation); + break; + } + case (dd::Y): { + dd_operation = localDD->multiply(localDD->makeGateDD(dd::Ymat, getNumberOfQubits(), target), dd_operation); + break; + } + case (dd::Z): { + dd_operation = localDD->multiply(localDD->makeGateDD(dd::Zmat, getNumberOfQubits(), target), dd_operation); + break; + } + default: { + throw std::runtime_error("Unknown noise operation '" + std::to_string(effect) + "'\n"); + } + } + } + } + return dd_operation; +} + +dd::Package::vEdge StochasticNoiseSimulator::MeasureOneCollapsingConcurrent(unsigned short index, + const std::unique_ptr& localDD, + dd::Package::vEdge local_root_edge, + std::mt19937_64& generator, + char* result, + bool assume_probability_normalization) { + //todo merge the function with MeasureOneCollapsing + std::map probsMone; + std::set visited_nodes2; + std::queue q; + + probsMone[local_root_edge.p] = dd::ComplexNumbers::mag2(local_root_edge.w); + visited_nodes2.insert(local_root_edge.p); + q.push(local_root_edge.p); + + while (q.front()->v != index) { + dd::Package::vNode* ptr = q.front(); + q.pop(); + double prob = probsMone[ptr]; + + if (!ptr->e.at(0).w.approximatelyZero()) { + const dd::fp tmp1 = prob * dd::ComplexNumbers::mag2(ptr->e.at(0).w); + + if (visited_nodes2.find(ptr->e.at(0).p) != visited_nodes2.end()) { + probsMone[ptr->e.at(0).p] = probsMone[ptr->e.at(0).p] + tmp1; + } else { + probsMone[ptr->e.at(0).p] = tmp1; + visited_nodes2.insert(ptr->e.at(0).p); + q.push(ptr->e.at(0).p); + } + } + + if (!ptr->e.at(1).w.approximatelyZero()) { + const dd::fp tmp1 = prob * dd::ComplexNumbers::mag2(ptr->e.at(1).w); + + if (visited_nodes2.find(ptr->e.at(1).p) != visited_nodes2.end()) { + probsMone[ptr->e.at(1).p] = probsMone[ptr->e.at(1).p] + tmp1; + } else { + probsMone[ptr->e.at(1).p] = tmp1; + visited_nodes2.insert(ptr->e.at(1).p); + q.push(ptr->e.at(1).p); + } + } + } + + dd::fp pzero{0}, pone{0}; + + if (assume_probability_normalization) { + while (!q.empty()) { + dd::Package::vNode* ptr = q.front(); + q.pop(); + + if (!ptr->e.at(0).w.approximatelyZero()) { + pzero += probsMone[ptr] * dd::ComplexNumbers::mag2(ptr->e.at(0).w); + } + + if (!ptr->e.at(1).w.approximatelyZero()) { + pone += probsMone[ptr] * dd::ComplexNumbers::mag2(ptr->e.at(1).w); + } + } + } else { + std::unordered_map probs; + assign_probs(root_edge, probs); + + while (!q.empty()) { + dd::Package::vNode* ptr = q.front(); + q.pop(); + + if (!ptr->e.at(0).w.approximatelyZero()) { + pzero += probsMone[ptr] * probs[ptr->e.at(0).p] * dd::ComplexNumbers::mag2(ptr->e.at(0).w); + } + + if (!ptr->e.at(1).w.approximatelyZero()) { + pone += probsMone[ptr] * probs[ptr->e.at(1).p] * dd::ComplexNumbers::mag2(ptr->e.at(1).w); + } + } + } + + if (std::abs(pzero + pone - 1) > epsilon) { + throw std::runtime_error("Numerical instability occurred during measurement: |alpha|^2 + |beta|^2 = " + std::to_string(pzero) + " + " + std::to_string(pone) + " = " + std::to_string(pzero + pone) + ", but should be 1!"); + } + + const dd::fp sum = pzero + pone; + + //std::pair probs = std::make_pair(pzero, pone); + dd::GateMatrix measure_m{ + dd::complex_zero, dd::complex_zero, + dd::complex_zero, dd::complex_zero}; + + std::uniform_real_distribution dist(0.0, 1.0L); + + dd::fp n = dist(generator); + dd::fp norm_factor; + + if (n < pzero / sum) { + measure_m[0] = dd::complex_one; + norm_factor = pzero; + *result = '0'; + } else { + measure_m[3] = dd::complex_one; + norm_factor = pone; + *result = '1'; + } + dd::Edge m_gate = localDD->makeGateDD(measure_m, getNumberOfQubits(), index); + + dd::Edge e = localDD->multiply(m_gate, local_root_edge); + + dd::Complex c = localDD->cn.getTemporary(std::sqrt(1 / norm_factor), 0); + dd::ComplexNumbers::mul(c, e.w, c); + e.w = localDD->cn.lookup(c); + localDD->incRef(e); + localDD->decRef(local_root_edge); + local_root_edge = e; + + return local_root_edge; +} + +dd::NoiseOperationKind StochasticNoiseSimulator::ReturnNoiseOperation(char i, double prob) const { + // if (forcedResult.size() > currentResult) { + // auto tmp = forcedResult[currentResult]; + // currentResult += 1; + // return tmp; + // } + switch (i) { + case 'D': { + if (prob >= 3 * noise_probability / 4) { + // printf("T"); + return dd::I; + } else if (prob < noise_probability / 4) { + // printf("X"); + return dd::X; + } else if (prob < noise_probability / 2) { + // printf("Y"); + return dd::Y; + } else { //if (prob < 3 * noise_probability / 4) { + // printf("Z"); + return dd::Z; + } + } + case 'A': { + return dd::AFalse; + } + case 'P': { + if (prob > noise_probability) { + // printf("B"); + return dd::I; + } else { + // printf("P"); + return dd::Z; // A phase flip is represented by a Z operation + } + } + default: + throw std::runtime_error(std::string{"Unknown noise effect '"} + i + "'"); + } +} + +void StochasticNoiseSimulator::setRecordedProperties(const std::string& input) { + std::string subString; + int list_begin = std::numeric_limits::min(); + int list_end = std::numeric_limits::min(); + for (char i: input) { + if (i == ' ') { + continue; + } + if (i == ',') { + if (list_begin > std::numeric_limits::min()) { + list_end = std::stoi(subString); + if (list_end > std::pow(2, getNumberOfQubits())) { + list_end = (int)std::pow(2, getNumberOfQubits()); + } + assert(list_begin < list_end); + for (int m = list_begin; m <= list_end; m++) { + recorded_properties.emplace_back(std::make_tuple(m, intToString(m))); + } + } else { + recorded_properties.emplace_back(std::make_tuple(std::stoi(subString), intToString(std::stoi(subString)))); + } + subString = ""; + list_begin = std::numeric_limits::min(); + list_end = std::numeric_limits::min(); + continue; + } + if (i == '-' && subString.length() > 0) { + list_begin = std::stoi(subString); + subString = ""; + continue; + } + subString += i; + } + if (list_begin > std::numeric_limits::min()) { + list_end = std::stoi(subString); + if (list_end > std::pow(2, getNumberOfQubits())) { + list_end = (int)std::pow(2, getNumberOfQubits()) - 1; + } + assert(list_begin < list_end); + for (int m = list_begin; m <= list_end; m++) { + recorded_properties.emplace_back(std::make_tuple(m, intToString(m))); + } + } else { + recorded_properties.emplace_back(std::make_tuple(std::stoi(subString), intToString(std::stoi(subString)))); + } +} + +std::string StochasticNoiseSimulator::intToString(long target_number) const { + if (target_number < 0) { + return "X"; + } + auto qubits = getNumberOfQubits(); + std::string path(qubits, '0'); + auto number = (unsigned long)target_number; + for (int i = 1; i <= qubits; i++) { + if (number % 2) { + path[qubits - i] = '1'; + } + number = number >> 1u; + } + return path; +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 339e679d..99431483 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -26,6 +26,7 @@ package_add_test(${PROJECT_NAME}_test ${CMAKE_CURRENT_SOURCE_DIR}/test_shor_sim.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_fast_shor_sim.cpp ${CMAKE_CURRENT_SOURCE_DIR}/test_grover_sim.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/test_stoch_noise_sim.cpp ) add_custom_command(TARGET ${PROJECT_NAME}_test diff --git a/test/test_stoch_noise_sim.cpp b/test/test_stoch_noise_sim.cpp new file mode 100644 index 00000000..8ad16e49 --- /dev/null +++ b/test/test_stoch_noise_sim.cpp @@ -0,0 +1,424 @@ +#include "StochasticNoiseSimulator.hpp" + +#include +#include + +/** + * These tests may have to be adjusted if something about the random-number generation changes. + */ + +TEST(StochNoiseSimTest, SingleOneQubitGateOnTwoQubitCircuit) { + auto quantumComputation = std::make_unique(2); + quantumComputation->emplace_back(2, 0, qc::X); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + + ASSERT_EQ(ddsim.getNumberOfOps(), 1); + + ddsim.Simulate(1); + + auto m = ddsim.MeasureAll(false); + + ASSERT_EQ("01", m); +} + +TEST(StochNoiseSimTest, ClassicControlledOp) { + auto quantumComputation = std::make_unique(2); + quantumComputation->emplace_back(2, 0, qc::X); + quantumComputation->emplace_back(2, 0, 0); + std::unique_ptr op(new qc::StandardOperation(2, 1, qc::X)); + auto classical_register = std::make_pair(0, 1); + quantumComputation->emplace_back(op, classical_register, 1); + + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + + auto m = ddsim.MeasureAll(false); + + ASSERT_EQ("11", m); +} + +TEST(StochNoiseSimTest, DestructiveMeasurementAll) { + auto quantumComputation = std::make_unique(2); + quantumComputation->emplace_back(2, 0, qc::H); + quantumComputation->emplace_back(2, 1, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + + const std::vector v_before = ddsim.getVector(); + ASSERT_EQ(v_before[0], v_before[1]); + ASSERT_EQ(v_before[0], v_before[2]); + ASSERT_EQ(v_before[0], v_before[3]); + + const std::string m = ddsim.MeasureAll(true); + const std::vector v_after = ddsim.getVector(); + const int i = std::stoi(m, nullptr, 2); + + ASSERT_EQ(v_after[i].r, 1.0); + ASSERT_EQ(v_after[i].i, 0.0); +} + +TEST(StochNoiseSimTest, DestructiveMeasurementOne) { + auto quantumComputation = std::make_unique(2); + quantumComputation->emplace_back(2, 0, qc::H); + quantumComputation->emplace_back(2, 1, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + + const char m = ddsim.MeasureOneCollapsing(0); + const std::vector v_after = ddsim.getVector(); + + if (m == '0') { + ASSERT_EQ(v_after[0], dd::complex_SQRT2_2); + ASSERT_EQ(v_after[2], dd::complex_SQRT2_2); + ASSERT_EQ(v_after[1], (dd::ComplexValue{0, 0})); + ASSERT_EQ(v_after[3], (dd::ComplexValue{0, 0})); + } else if (m == '1') { + ASSERT_EQ(v_after[0], (dd::ComplexValue{0, 0})); + ASSERT_EQ(v_after[2], (dd::ComplexValue{0, 0})); + ASSERT_EQ(v_after[1], dd::complex_SQRT2_2); + ASSERT_EQ(v_after[3], dd::complex_SQRT2_2); + } else { + FAIL() << "Measurement result not in {0,1}!"; + } +} + +TEST(StochNoiseSimTest, DestructiveMeasurementOneArbitraryNormalization) { + auto quantumComputation = std::make_unique(2); + quantumComputation->emplace_back(2, 0, qc::H); + quantumComputation->emplace_back(2, 1, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + + char m = 'X'; + std::mt19937_64 gen{}; + + ddsim.root_edge = ddsim.MeasureOneCollapsingConcurrent(0, ddsim.dd, ddsim.root_edge, gen, &m, false); + + const std::vector v_after = ddsim.getVector(); + + for (auto const& e: v_after) { + std::cout << e << " "; + } + std::cout << "\n"; + + if (m == '0') { + ASSERT_EQ(v_after[0], dd::complex_SQRT2_2); + ASSERT_EQ(v_after[2], dd::complex_SQRT2_2); + ASSERT_EQ(v_after[1], (dd::ComplexValue{0, 0})); + ASSERT_EQ(v_after[3], (dd::ComplexValue{0, 0})); + } else if (m == '1') { + ASSERT_EQ(v_after[0], (dd::ComplexValue{0, 0})); + ASSERT_EQ(v_after[2], (dd::ComplexValue{0, 0})); + ASSERT_EQ(v_after[1], dd::complex_SQRT2_2); + ASSERT_EQ(v_after[3], dd::complex_SQRT2_2); + } else { + FAIL() << "Measurement result not in {0,1}!"; + } +} + +TEST(StochNoiseSimTest, ApproximateByFidelity) { + auto quantumComputation = std::make_unique(3); + quantumComputation->emplace_back(3, 0, qc::H); + quantumComputation->emplace_back(3, 1, qc::H); + quantumComputation->emplace_back(3, dd::Controls{dd::Control{0}, dd::Control{1}}, 2, qc::X); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + + ASSERT_EQ(ddsim.getActiveNodeCount(), 6); + + double resulting_fidelity = ddsim.ApproximateByFidelity(0.3, false, true); + + ASSERT_EQ(ddsim.getActiveNodeCount(), 4); + ASSERT_DOUBLE_EQ(resulting_fidelity, 0.75); //equal up to 4 ULP +} + +TEST(StochNoiseSimTest, ApproximateBySampling) { + auto quantumComputation = std::make_unique(3); + quantumComputation->emplace_back(3, 0, qc::H); + quantumComputation->emplace_back(3, 1, qc::H); + quantumComputation->emplace_back(3, dd::Controls{dd::Control{0}, dd::Control{1}}, 2, qc::X); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + + ASSERT_EQ(ddsim.getActiveNodeCount(), 6); + + double resulting_fidelity = ddsim.ApproximateBySampling(1, 0, true); + + ASSERT_EQ(ddsim.getActiveNodeCount(), 3); + ASSERT_LE(resulting_fidelity, 0.75); // the least contributing path has .25 +} + +TEST(StochNoiseSimTest, Reordering) { + auto quantumComputation = std::make_unique(3); + quantumComputation->emplace_back(3, 0, qc::H); + quantumComputation->emplace_back(3, 1, qc::H); + quantumComputation->emplace_back(3, std::vector{0, 1, 2}, qc::OpType::Barrier); + quantumComputation->emplace_back(3, dd::Controls{dd::Control{0}, dd::Control{1}}, 2, qc::X); + + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); +} + +TEST(StochNoiseSimTest, SimulateAdder4TrackFidelityWithNoise) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, std::string("APD"), 0.01, 30000, 1, 1, "-3-1000"); + auto m = ddsim.StochSimulate(); + + std::cout << ddsim.getName() << "\n"; + + EXPECT_GE(m.find("1000")->second, 0.14); + EXPECT_LE(m.find("1000")->second, 0.16); + EXPECT_GE(m.find("1001")->second, 0.54); + EXPECT_LE(m.find("1001")->second, 0.56); + EXPECT_GE(std::stod(ddsim.AdditionalStatistics().at("final_fidelity")), 0.54); + EXPECT_LE(std::stod(ddsim.AdditionalStatistics().at("final_fidelity")), 0.56); +} + +TEST(StochNoiseSimTest, SimulateClassicControlledOpWithError) { + auto quantumComputation = std::make_unique(2); + quantumComputation->emplace_back(2, 0, qc::X); + quantumComputation->emplace_back(2, 0, 0); + quantumComputation->emplace_back(2, 0, qc::H); + std::unique_ptr op(new qc::StandardOperation(2, 1, qc::X)); + auto classical_register = std::make_pair(0, 1); + quantumComputation->emplace_back(op, classical_register, 1); + + StochasticNoiseSimulator ddsim(quantumComputation, std::string("APD"), 0.02, 30000, 1, 1, "-1-1000"); + auto m = ddsim.StochSimulate(); + + EXPECT_GE(m.find("00")->second, 0.03); + EXPECT_LE(m.find("00")->second, 0.06); + EXPECT_GE(m.find("01")->second, 0.46); + EXPECT_LE(m.find("01")->second, 0.48); + EXPECT_GE(m.find("10")->second, 0.03); + EXPECT_LE(m.find("10")->second, 0.05); + EXPECT_GE(m.find("11")->second, 0.42); + EXPECT_LE(m.find("11")->second, 0.44); +} + +TEST(StochNoiseSimTest, SimulateAdder4WithoutNoise) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, 1, 1); + ddsim.Simulate(1); + auto m = ddsim.MeasureAll(false); + ASSERT_EQ("1001", m); +} + +TEST(StochNoiseSimTest, SimulateAdder4WithDecoherenceAndGateError) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, std::string("APD"), 0.1, 30000, 1, 1, "-1-1000"); + auto m = ddsim.StochSimulate(); + + EXPECT_GE(m.find("0000")->second, 0.19); + EXPECT_LE(m.find("0000")->second, 0.23); + EXPECT_GE(m.find("0001")->second, 0.13); + EXPECT_LE(m.find("0001")->second, 0.15); + EXPECT_GE(m.find("0100")->second, 0.08); + EXPECT_LE(m.find("0100")->second, 0.10); +} + +TEST(StochNoiseSimTest, SimulateAdder4WithDecoherenceAndGateErrorSelectedProperties) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, std::vector{0, 1, 2}, qc::OpType::Barrier); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, std::string("APD"), 0.1, 30000, 1, 1, " -3-500,501,502"); + auto m = ddsim.StochSimulate(); + + EXPECT_GE(m.find("0000")->second, 0.19); + EXPECT_LE(m.find("0000")->second, 0.23); + EXPECT_GE(m.find("0001")->second, 0.13); + EXPECT_LE(m.find("0001")->second, 0.15); + EXPECT_GE(m.find("0100")->second, 0.08); + EXPECT_LE(m.find("0100")->second, 0.10); +} + +TEST(StochNoiseSimTest, SimulateAdder4WithDecoherenceError) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, std::string("AP"), 0.3, 30000, 1, 1, "0-1000"); + auto m = ddsim.StochSimulate(); + + EXPECT_GE(m.find("0000")->second, 0.78); + EXPECT_LE(m.find("0000")->second, 0.80); + EXPECT_GE(m.find("0001")->second, 0.18); + EXPECT_LE(m.find("0001")->second, 0.20); +} + +TEST(StochNoiseSimTest, SimulateAdder4WithDepolarizationError) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, std::string("D"), 0.01, 30000, 1, 1, "0-1000"); + auto m = ddsim.StochSimulate(); + + EXPECT_GE(m.find("1001")->second, 0.81); + EXPECT_LE(m.find("1001")->second, 0.84); + EXPECT_GE(m.find("0000")->second, 0.01); + EXPECT_LE(m.find("0000")->second, 0.03); +} + +TEST(StochNoiseSimTest, SimulateAdder4WithNoiseAndApproximation) { + auto quantumComputation = std::make_unique(4); + quantumComputation->emplace_back(4, 0, qc::X); + quantumComputation->emplace_back(4, 1, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::T); + quantumComputation->emplace_back(4, 1, qc::T); + quantumComputation->emplace_back(4, 2, qc::T); + quantumComputation->emplace_back(4, 3, qc::Tdag); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{1}}, 2, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 0, qc::Tdag); + quantumComputation->emplace_back(4, 1, qc::Tdag); + quantumComputation->emplace_back(4, 2, qc::Tdag); + quantumComputation->emplace_back(4, 3, qc::T); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{0}}, 1, qc::X); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{2}}, 3, qc::X); + quantumComputation->emplace_back(4, 3, qc::S); + quantumComputation->emplace_back(4, dd::Controls{dd::Control{3}}, 0, qc::X); + quantumComputation->emplace_back(4, 3, qc::H); + StochasticNoiseSimulator ddsim(quantumComputation, std::string("APD"), 0.01, 30000, 1, 0.9, "-3-1000"); + auto m = ddsim.StochSimulate(); + + EXPECT_DOUBLE_EQ(std::stod(ddsim.AdditionalStatistics().at("step_fidelity")), 0.9); + EXPECT_GT(std::stoi(ddsim.AdditionalStatistics().at("approximation_runs")), 0); +} \ No newline at end of file