From 8616db6698b297a69c7f422acfe2266223d28d1a Mon Sep 17 00:00:00 2001 From: AlexandreF-1qbit <76115575+AlexandreF-1qbit@users.noreply.github.com> Date: Wed, 28 Sep 2022 13:46:28 -0400 Subject: [PATCH] Histogram class (#217) * Histogram class, post-selection, resampling and aggregation methods --- examples/overview_endtoend.ipynb | 322 ++++++++---------- tangelo/toolboxes/post_processing/__init__.py | 1 + .../toolboxes/post_processing/histogram.py | 227 ++++++++++++ .../post_processing/tests/test_histogram.py | 200 +++++++++++ 4 files changed, 574 insertions(+), 176 deletions(-) create mode 100644 tangelo/toolboxes/post_processing/histogram.py create mode 100644 tangelo/toolboxes/post_processing/tests/test_histogram.py diff --git a/examples/overview_endtoend.ipynb b/examples/overview_endtoend.ipynb index e4a329b1a..d42da6b93 100644 --- a/examples/overview_endtoend.ipynb +++ b/examples/overview_endtoend.ipynb @@ -53,10 +53,8 @@ }, { "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, + "execution_count": 1, + "metadata": {}, "outputs": [], "source": [ "# Pretty printer for more readable outputs\n", @@ -101,9 +99,7 @@ { "cell_type": "code", "execution_count": 2, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -112,31 +108,35 @@ "{'active_occupied': [0, 1, 2, 3, 4],\n", " 'active_virtual': [5, 6, 7, 8, 9],\n", " 'basis': 'minao',\n", + " 'ecp': {},\n", " 'frozen_occupied': [],\n", " 'frozen_orbitals': None,\n", " 'frozen_virtual': [],\n", - " 'mean_field': ,\n", - " 'mf_energy': -5.27727228441684,\n", + " 'mean_field': RHF object of ,\n", + " 'mf_energy': -5.277272284416819,\n", " 'mo_energies': array([-0.71250664, -0.63494046, -0.63486503, -0.4162918 , -0.41611742,\n", " 0.07821244, 0.07835203, 0.50868475, 0.50889667, 0.75894782]),\n", " 'mo_occ': array([2., 2., 2., 2., 2., 0., 0., 0., 0., 0.]),\n", + " 'mo_symm_ids': None,\n", + " 'mo_symm_labels': None,\n", " 'n_atoms': 10,\n", " 'n_electrons': 10,\n", - " 'n_min_orbitals': 10,\n", + " 'n_min_orbitals': 40,\n", " 'n_mos': 10,\n", " 'n_sos': 20,\n", " 'q': 0,\n", " 'spin': 0,\n", - " 'xyz': [['H', (0.0, 1.78, 0.0)],\n", - " ['H', (-1.046, 1.44, 0.0)],\n", - " ['H', (-1.693, 0.55, 0.0)],\n", - " ['H', (-1.693, -0.55, 0.0)],\n", - " ['H', (-1.046, -1.44, 0.0)],\n", - " ['H', (0.0, -1.78, 0.0)],\n", - " ['H', (1.046, -1.44, 0.0)],\n", - " ['H', (1.693, -0.55, 0.0)],\n", - " ['H', (1.693, 0.55, 0.0)],\n", - " ['H', (1.046, 1.44, 0.0)]]}\n" + " 'symmetry': False,\n", + " 'xyz': [('H', (0.0, 1.78, 0.0)),\n", + " ('H', (-1.046, 1.44, 0.0)),\n", + " ('H', (-1.6929999999999996, 0.55, 0.0)),\n", + " ('H', (-1.6929999999999996, -0.55, 0.0)),\n", + " ('H', (-1.046, -1.44, 0.0)),\n", + " ('H', (0.0, -1.78, 0.0)),\n", + " ('H', (1.046, -1.44, 0.0)),\n", + " ('H', (1.6929999999999996, -0.55, 0.0)),\n", + " ('H', (1.6929999999999996, 0.55, 0.0)),\n", + " ('H', (1.046, 1.44, 0.0))]}\n" ] } ], @@ -174,9 +174,7 @@ { "cell_type": "code", "execution_count": 3, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -216,9 +214,7 @@ { "cell_type": "code", "execution_count": 4, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -271,9 +267,7 @@ { "cell_type": "code", "execution_count": 5, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -288,7 +282,7 @@ " 'vqe_variational_parameters': 14}\n", "\n", "DMET-VQE-UCCSD, 10 fragments \n", - "{'circuit_2qubit_gates': 4, 'circuit_gates': 20, 'circuit_var_gates': 4, 'circuit_width': 2, 'qubit_hamiltonian_terms': 9, 'vqe_variational_parameters': 2}\n", + "{'circuit_2qubit_gates': 4, 'circuit_gates': 22, 'circuit_var_gates': 4, 'circuit_width': 2, 'qubit_hamiltonian_terms': 9, 'vqe_variational_parameters': 2}\n", "\n" ] } @@ -348,9 +342,7 @@ { "cell_type": "code", "execution_count": 6, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "dmet_energy = dmet_solver.simulate()" @@ -359,16 +351,14 @@ { "cell_type": "code", "execution_count": 7, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - " DMET-VQE 10 fragments: -5.40145\n", - " Difference with FCI: 8.625559E-03\n" + " DMET-VQE 10 fragments: -5.40161\n", + " Difference with FCI: 8.464849E-03\n" ] } ], @@ -392,36 +382,36 @@ { "cell_type": "code", "execution_count": 8, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Circuit object. Size 20 \n", + "Circuit object. Size 22 \n", "\n", - "RX target : 0 parameter : 1.5707963267948966\n", - "RZ target : 0 parameter : 0.0005896947252139878\t (variational)\n", - "RX target : 0 parameter : 10.995574287564276\n", - "RX target : 1 parameter : 1.5707963267948966\n", - "RZ target : 1 parameter : 0.0005896947252139878\t (variational)\n", - "RX target : 1 parameter : 10.995574287564276\n", - "RX target : 0 parameter : 1.5707963267948966\n", - "H target : 1 \n", - "CNOT target : 1 control : 0 \n", - "RZ target : 1 parameter : 1.6568407597743045\t (variational)\n", - "CNOT target : 1 control : 0 \n", - "H target : 1 \n", - "RX target : 0 parameter : 10.995574287564276\n", - "H target : 0 \n", - "RX target : 1 parameter : 1.5707963267948966\n", - "CNOT target : 1 control : 0 \n", - "RZ target : 1 parameter : 1.6568407597743045\t (variational)\n", - "CNOT target : 1 control : 0 \n", - "RX target : 1 parameter : 10.995574287564276\n", - "H target : 0 \n", + "X target : [0] \n", + "X target : [1] \n", + "RX target : [0] parameter : 1.5707963267948966\n", + "RZ target : [0] parameter : 3.1415958073389\t (variational)\n", + "RX target : [0] parameter : -1.5707963267948966\n", + "RX target : [1] parameter : 1.5707963267948966\n", + "RZ target : [1] parameter : 3.1415958073389\t (variational)\n", + "RX target : [1] parameter : -1.5707963267948966\n", + "RX target : [0] parameter : 1.5707963267948966\n", + "H target : [1] \n", + "CNOT target : [1] control : [0] \n", + "RZ target : [1] parameter : 1.658627529912278\t (variational)\n", + "CNOT target : [1] control : [0] \n", + "H target : [1] \n", + "RX target : [0] parameter : -1.5707963267948966\n", + "H target : [0] \n", + "RX target : [1] parameter : 1.5707963267948966\n", + "CNOT target : [1] control : [0] \n", + "RZ target : [1] parameter : 1.658627529912278\t (variational)\n", + "CNOT target : [1] control : [0] \n", + "RX target : [1] parameter : -1.5707963267948966\n", + "H target : [0] \n", "\n" ] } @@ -449,9 +439,7 @@ { "cell_type": "code", "execution_count": 9, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -512,9 +500,7 @@ { "cell_type": "code", "execution_count": 10, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -552,38 +538,38 @@ { "cell_type": "code", "execution_count": 11, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Circuit object. Size 22 \n", + "Circuit object. Size 24 \n", "\n", - "RX target : 0 parameter : 1.5707963267948966\n", - "RZ target : 0 parameter : 0.0005896947252139878\t (variational)\n", - "RX target : 0 parameter : 10.995574287564276\n", - "RX target : 1 parameter : 1.5707963267948966\n", - "RZ target : 1 parameter : 0.0005896947252139878\t (variational)\n", - "RX target : 1 parameter : 10.995574287564276\n", - "RX target : 0 parameter : 1.5707963267948966\n", - "H target : 1 \n", - "CNOT target : 1 control : 0 \n", - "RZ target : 1 parameter : 1.6568407597743045\t (variational)\n", - "CNOT target : 1 control : 0 \n", - "H target : 1 \n", - "RX target : 0 parameter : 10.995574287564276\n", - "H target : 0 \n", - "RX target : 1 parameter : 1.5707963267948966\n", - "CNOT target : 1 control : 0 \n", - "RZ target : 1 parameter : 1.6568407597743045\t (variational)\n", - "CNOT target : 1 control : 0 \n", - "RX target : 1 parameter : 10.995574287564276\n", - "H target : 0 \n", - "RX target : 0 parameter : 1.5707963267948966\n", - "RX target : 1 parameter : 1.5707963267948966\n", + "X target : [0] \n", + "X target : [1] \n", + "RX target : [0] parameter : 1.5707963267948966\n", + "RZ target : [0] parameter : 3.1415958073389\t (variational)\n", + "RX target : [0] parameter : -1.5707963267948966\n", + "RX target : [1] parameter : 1.5707963267948966\n", + "RZ target : [1] parameter : 3.1415958073389\t (variational)\n", + "RX target : [1] parameter : -1.5707963267948966\n", + "RX target : [0] parameter : 1.5707963267948966\n", + "H target : [1] \n", + "CNOT target : [1] control : [0] \n", + "RZ target : [1] parameter : 1.658627529912278\t (variational)\n", + "CNOT target : [1] control : [0] \n", + "H target : [1] \n", + "RX target : [0] parameter : -1.5707963267948966\n", + "H target : [0] \n", + "RX target : [1] parameter : 1.5707963267948966\n", + "CNOT target : [1] control : [0] \n", + "RZ target : [1] parameter : 1.658627529912278\t (variational)\n", + "CNOT target : [1] control : [0] \n", + "RX target : [1] parameter : -1.5707963267948966\n", + "H target : [0] \n", + "RX target : [0] parameter : 1.5707963267948966\n", + "RX target : [1] parameter : 1.5707963267948966\n", "\n" ] } @@ -619,9 +605,7 @@ { "cell_type": "code", "execution_count": 12, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -656,9 +640,7 @@ { "cell_type": "code", "execution_count": 13, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "n_shots = 10000" @@ -704,9 +686,23 @@ "from tangelo.linq.qpu_connection import QEMISTCloudConnection\n", "\n", "qcloud_connection = QEMISTCloudConnection()\n", - "\n", - "circuit_YY = quantum_circuit[((0, \"Y\"), (1, \"Y\"))]\n", - "\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "circuit_YY = quantum_circuit[((0, \"Y\"), (1, \"Y\"))]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```python\n", "price_estimates = qcloud_connection.job_estimate(circuit_YY, n_shots=n_shots)\n", "print(price_estimates)\n", "\n", @@ -714,13 +710,13 @@ "\n", "# This two commands would respecfully submit the job to the quantum device and make a blocking call\n", "# to retrieve the results, through a job ID returned by QEMIST Cloud\n", - "# job_id = qcloud_connection.job_submit(circuit_YY, n_shots=n_shots, backend=backend)\n", - "# freqs, raw_data = qcloud_connection.job_results(job_id)\n", + "job_id = qcloud_connection.job_submit(circuit_YY, n_shots=n_shots, backend=backend)\n", + "freqs, raw_data = qcloud_connection.job_results(job_id)\n", "```\n", "\n", "Output:\n", "```output\n", - "{'braket_ionq': 100.3, 'braket_rigetti': 3.8}\n", + "{'arn:aws:braket:::device/qpu/ionq/ionQdevice': 100.3, 'arn:aws:braket:::device/quantum-simulator/amazon/sv1': 3.8, 'arn:aws:braket:us-west-1::device/qpu/rigetti/Aspen-M-1': 3.8}\n", "```" ] }, @@ -738,21 +734,19 @@ { "cell_type": "code", "execution_count": 15, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "T : | 0 | 1 | 2 | 3 |4| 5 |6| 7 | 8 |9| 10 |11| 12 | 13 |\n", - " \n", - "q0 : -Rx(1.57)-Rz(0.00059)-Rx(11)-Rx(1.57)-C----------C-Rx(11)-H--------C----------C--H------Rx(1.57)-\n", - " | | | | \n", - "q1 : -Rx(1.57)-Rz(0.00059)-Rx(11)-H--------X-Rz(1.66)-X-H------Rx(1.57)-X-Rz(1.66)-X--Rx(11)-Rx(1.57)-\n", + "T : |0| 1 | 2 | 3 | 4 |5| 6 |7| 8 | 9 |10| 11 |12| 13 | 14 |\n", + " \n", + "q0 : -X-Rx(1.57)-Rz(3.14)-Rx(-1.57)-Rx(1.57)-C----------C-Rx(-1.57)-H--------C-----------C--H---------Rx(1.57)-\n", + " | | | | \n", + "q1 : -X-Rx(1.57)-Rz(3.14)-Rx(-1.57)-H--------X-Rz(1.66)-X-H---------Rx(1.57)-X--Rz(1.66)-X--Rx(-1.57)-Rx(1.57)-\n", "\n", - "T : | 0 | 1 | 2 | 3 |4| 5 |6| 7 | 8 |9| 10 |11| 12 | 13 |\n" + "T : |0| 1 | 2 | 3 | 4 |5| 6 |7| 8 | 9 |10| 11 |12| 13 | 14 |\n" ] } ], @@ -773,17 +767,15 @@ { "cell_type": "code", "execution_count": 16, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "0: ───I───Rx(0.5π)───Rz(0)───Rx(-0.5π)───Rx(0.5π)───@────────────────@───Rx(-0.5π)───H──────────@────────────────@───H───────────Rx(0.5π)───\n", - " │ │ │ │\n", - "1: ───I───Rx(0.5π)───Rz(0)───Rx(-0.5π)───H──────────X───Rz(0.527π)───X───H───────────Rx(0.5π)───X───Rz(0.527π)───X───Rx(-0.5π)───Rx(0.5π)───\n" + "0: ───I───X───Rx(0.5π)───Rz(π)───Rx(-0.5π)───Rx(0.5π)───@────────────────@───Rx(-0.5π)───H──────────@────────────────@───H───────────Rx(0.5π)───\n", + " │ │ │ │\n", + "1: ───I───X───Rx(0.5π)───Rz(π)───Rx(-0.5π)───H──────────X───Rz(0.528π)───X───H───────────Rx(0.5π)───X───Rz(0.528π)───X───Rx(-0.5π)───Rx(0.5π)───\n" ] } ], @@ -810,9 +802,7 @@ { "cell_type": "code", "execution_count": 17, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "from tangelo.linq import Simulator\n", @@ -829,15 +819,13 @@ { "cell_type": "code", "execution_count": 18, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "{'00': 0.2879, '01': 0.2126, '10': 0.2076, '11': 0.2919}\n" + "{'00': 0.2905, '01': 0.2082, '10': 0.2062, '11': 0.2951}\n" ] } ], @@ -876,9 +864,7 @@ { "cell_type": "code", "execution_count": 19, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "# DMET-QCC experimental frequencies.\n", @@ -891,9 +877,7 @@ { "cell_type": "code", "execution_count": 20, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "# Derived thanks to observations after compilation and optimization: similar to XZ\n", @@ -904,21 +888,21 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can then compute the expectation values of the different qubit operators of interest by using these frequencies. For simplicity, we manually compute all 9 of them, from the 5 histograms. In our case, since all circuits have been run with the same number of shots, it is pretty straightforward. Note that because the term grouping based on qubit-wise commutativity is not unique, several of these circuits could be used to compute the output of a different computational basis: in our case, it means that we overall have 20,000 shots worth of data by combining the two relevant histograms for certain entries:" + "We can then compute the expectation values of the different qubit operators of interest by using these frequencies. For simplicity, we manually compute all 9 of them, from the 5 histograms. Note that because the term grouping based on qubit-wise commutativity is not unique, several of these circuits could be used to compute the output of a different computational basis: in our case, it means that we overall have 20,000 shots worth of data by aggregating the two relevant histograms for certain entries,\n", + "\n", + "We provide a `Histogram` helper class and the `aggregate_histogram` function to facilitate this kind of operation. This data-structure encapsulates the dictionary of frequencies and the number of shots used to build it, which is used as a \\\"weight\\\" when aggregating with other histograms. If we keep track of the number of shots used for each computational basis (with a dictionary for example) then it is pretty straightforward to leverage that helper class. In our case it is even easier: all jobs used the name number of shots." ] }, { "cell_type": "code", "execution_count": 21, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "{((0, 'X'),): 0.0012000000000000344,\n", + "{((0, 'X'),): 0.0011999999999999789,\n", " ((0, 'X'), (1, 'X')): -0.17740000000000003,\n", " ((0, 'X'), (1, 'Z')): -0.0012000000000000344,\n", " ((0, 'Y'), (1, 'Y')): 0.17379999999999998,\n", @@ -931,21 +915,18 @@ } ], "source": [ - "def get_average_histogram(hists):\n", - " \"\"\" Compute the average of histograms provided as a list, assume identical weights (n_shots) for all of them \"\"\"\n", - " from collections import Counter\n", - " \n", - " sum_hist = sum([Counter(hist) for hist in hists], start=Counter())\n", - " avg_hist = {k:v/len(hists) for k,v in sum_hist.items()}\n", - " return avg_hist\n", + "from tangelo.toolboxes.post_processing import Histogram, aggregate_histograms\n", "\n", - "freq_dict[(0, 'Z'),] = get_average_histogram([freq_dict[(0, 'Z'), (1, 'Z')], freq_dict[(0, 'Z'), (1, 'X')]])\n", - "freq_dict[(1, 'Z'),] = get_average_histogram([freq_dict[(0, 'Z'), (1, 'Z')], freq_dict[(0, 'X'), (1, 'Z')]])\n", - "freq_dict[(0, 'X'),] = get_average_histogram([freq_dict[(0, 'X'), (1, 'Z')], freq_dict[(0, 'X'), (1, 'X')]])\n", - "freq_dict[(1, 'X'),] = get_average_histogram([freq_dict[(0, 'Z'), (1, 'X')], freq_dict[(0, 'X'), (1, 'X')]])\n", + "# Turn dictionary of frequencies into Histogram objects\n", + "freq_hists = {k: Histogram(freqs, n_shots) for k, freqs in freq_dict.items()}\n", "\n", - "expectation_values = {term: Simulator.get_expectation_value_from_frequencies_oneterm(term, hist) \n", - " for term, hist in freq_dict.items()}\n", + "# Compute the other Histogram objects by aggregating the ones obtained in the experiment\n", + "freq_hists[((0, 'Z'),)] = aggregate_histograms(freq_hists[((0, 'Z'), (1, 'Z'))], freq_hists[((0, 'Z'), (1, 'X'))])\n", + "freq_hists[((1, 'Z'),)] = aggregate_histograms(freq_hists[((0, 'Z'), (1, 'Z'))], freq_hists[((0, 'X'), (1, 'Z'))])\n", + "freq_hists[((0, 'X'),)] = aggregate_histograms(freq_hists[((0, 'X'), (1, 'Z'))], freq_hists[((0, 'X'), (1, 'X'))])\n", + "freq_hists[((1, 'X'),)] = aggregate_histograms(freq_hists[((0, 'Z'), (1, 'X'))], freq_hists[((0, 'X'), (1, 'X'))])\n", + "\n", + "expectation_values = {term: hist.get_expectation_value(term) for term, hist in freq_hists.items()}\n", "pprint(expectation_values)" ] }, @@ -966,9 +947,7 @@ { "cell_type": "code", "execution_count": 22, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -1031,16 +1010,14 @@ { "cell_type": "code", "execution_count": 23, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - " DMET-VQE QCC 1 fragment: -0.53827\n", - " Difference with FCI: 2.736958E-03\n" + " DMET-VQE QCC 1 fragment: -0.53820\n", + " Difference with FCI: 2.808912E-03\n" ] } ], @@ -1101,25 +1078,23 @@ { "cell_type": "code", "execution_count": 24, - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - " DMET-VQE QCC 1 fragment: -0.53893\n", - " Difference with FCI: 2.081728E-03\n" + " DMET-VQE QCC 1 fragment: -0.53885\n", + " Difference with FCI: 2.153904E-03\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "/home/valentin/Desktop/virtualenvs/tangelo_oct21/lib/python3.8/site-packages/tangelo/toolboxes/post_processing/mc_weeny_rdm_purification.py:68: ComplexWarning: Casting complex values to real discards the imaginary part\n", + "/home/alex/Codes/Tangelo/tangelo/toolboxes/post_processing/mc_weeny_rdm_purification.py:68: ComplexWarning: Casting complex values to real discards the imaginary part\n", " rdm1_np_temp[i, j] += D_matrix2_final[i, j, k, k]\n", - "/home/valentin/Desktop/virtualenvs/tangelo_oct21/lib/python3.8/site-packages/tangelo/toolboxes/post_processing/mc_weeny_rdm_purification.py:74: ComplexWarning: Casting complex values to real discards the imaginary part\n", + "/home/alex/Codes/Tangelo/tangelo/toolboxes/post_processing/mc_weeny_rdm_purification.py:74: ComplexWarning: Casting complex values to real discards the imaginary part\n", " rdm2_np[i//2, j//2, k//2, l//2] += D_matrix2_final[i, j, k, l]\n" ] } @@ -1140,15 +1115,13 @@ "\n", "Experimental data requires a measure of its uncertainty. As it is often prohibitively expensive to collect large amount of data on quantum computers for the purpose of estimating uncertainty, we generate statistics from our dataset using an established method called bootstrapping [[Efron, B. et al., 1994](https://books.google.ca/books/about/An_Introduction_to_the_Bootstrap.html?id=gLlpIUxRntoC&redir_esc=y)]. For each histogram obtained from our experiment, we resample with replacement from that distribution to generate new histograms of the same sample size. We then use these histograms to calculate a new set of expectation values, RDMs, and total energies. This process is repeated many times, and from those outcomes we calculate the average energy and standard deviation of our experiment.\n", "\n", - "In our publication, this process was repeated 10,000 times and led to a result of 0.540 ± 0.007 Hartree." + "In our publication, this process was repeated 10,000 times and led to a result of -0.540 ±0.007 Hartree." ] }, { "cell_type": "code", "execution_count": 25, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -1159,19 +1132,16 @@ } ], "source": [ - "from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies\n", - "\n", "# Bootstrap method.\n", "fragment_energies = list()\n", "\n", "for n in range(10000): # Was 10000 in our publication\n", " \n", " # Step 1-2: draw random bootstrap sample and construct new histograms.\n", - " resample_freq = {term: get_resampled_frequencies(freq, n_shots) for term, freq in freq_dict.items()}\n", + " resample_freq = {term: hist.resample(n_shots) for term, hist in freq_hists.items()}\n", " \n", " # Step 3-4: compute expectation values.\n", - " expectation_values = {term: Simulator.get_expectation_value_from_frequencies_oneterm(term, hist) \n", - " for term, hist in resample_freq.items()}\n", + " expectation_values = {term: hist.get_expectation_value(term) for term, hist in resample_freq.items()}\n", " \n", " # Step 5: construct 1- and 2-RDMs.\n", " onerdm, twordm, _, _ = compute_rdms(fragment, expectation_values)\n", diff --git a/tangelo/toolboxes/post_processing/__init__.py b/tangelo/toolboxes/post_processing/__init__.py index dab23877f..5e801a8fd 100644 --- a/tangelo/toolboxes/post_processing/__init__.py +++ b/tangelo/toolboxes/post_processing/__init__.py @@ -12,5 +12,6 @@ # See the License for the specific language governing permissions and # limitations under the License. +from .histogram import Histogram, aggregate_histograms, filter_hist from .mc_weeny_rdm_purification import mcweeny_purify_2rdm from .extrapolation import diis, richardson diff --git a/tangelo/toolboxes/post_processing/histogram.py b/tangelo/toolboxes/post_processing/histogram.py new file mode 100644 index 000000000..8c4c32a5b --- /dev/null +++ b/tangelo/toolboxes/post_processing/histogram.py @@ -0,0 +1,227 @@ +# Copyright 2021 Good Chemistry Company. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" This module provides a Histogram class and functions to manipulate them, +in order to facilitate post-processing of quantum experiments. +""" + +from collections import Counter + +from tangelo.linq import Simulator +from tangelo.toolboxes.post_processing.bootstrapping import get_resampled_frequencies + + +class Histogram: + """Class to provide useful tools helping redundant tasks when analyzing + data from an experiment. The expected data input is an histogram of + bitstrings ("010..."), where 0 (resp. 1) correspond to the |0> (resp. |1>) + measured state in the computational basis. The outcomes can refer to either + the shot counts for each computational basis, or the corresponding + probability. + + The internal representation is kept as an histogram of shots. Normalized + quantities are accessible via the Histogram "frequencies" properties. + Bitstrings are stored with the lowest-significant qubit first (lsq_first) + ordering. + + Args: + outcomes (dict of string: int or float): Results in the format of bitstring: + outcome, where outcome can be a number of shots a probability. + n_shots (int): Self-explanatory. If it is greater than 0, the class + considers that probabilities are provided. + msq_first (bool): Bit ordering. For example, 011 (msq_first) = 110 + (lsq_first). This depends on the hardware provider convention. + epsilon (float): When probabilities are provided, this parameter is used + to check if the sum is within [1-epsilon, 1+epsilon]. Default is + 1e-2. + """ + + def __init__(self, outcomes, n_shots=0, msq_first=False, epsilon=1e-2): + + # Error: bitstrings not of consistent length. + lengths_bitstrings = set([len(k) for k in outcomes.keys()]) + if len(lengths_bitstrings) > 1: + raise ValueError(f"Entries in outcomes dictionary not of consistent length. Please check your input.") + + if n_shots > 0: + # Flags histograms whose frequencies do not add to 1. + sum_values = sum(outcomes.values()) + if abs(sum_values-1) > epsilon: + raise ValueError(f"Sum of Histogram frequencies ({sum_values}) differ from 1. by more than an epsilon ({epsilon})." + "Please adjust the value of epsilon or adjust shot data.") + + outcomes = {k: round(v*n_shots) for k, v in outcomes.items()} + elif n_shots < 0: + raise ValueError(f"The number of shots provided ({n_shots}) must be a positive value.") + + self.counts = outcomes.copy() + + if msq_first: + self.counts = {k[::-1]: v for k, v in self.counts.items()} + + def __repr__(self): + """Output the string representation.""" + return f"{self.counts}" + + def __eq__(self, other): + """Two histograms are equal if they have same counts, i.e. same + bitstrings and outcomes. + """ + return (self.counts == other.counts) + + def __add__(self, other): + """The aggregate_histograms function is used to add two Histogram objects.""" + return aggregate_histograms(self, other) + + def __iadd__(self, other): + new_histogram = self + other + self.__dict__ = new_histogram.__dict__ + return self + + @property + def n_shots(self): + """Return the number of shots encoded in the Histogram. + + Returns: + int: Self-explanatory. + """ + return sum(self.counts.values()) + + @property + def n_qubits(self): + """The length of the bitstrings represents the number of qubits." + + Returns: + int: Self-explanatory. + """ + return len(next(iter(self.counts))) + + @property + def frequencies(self): + """The frequencies are normalized counts vs the number of shots. + + Returns: + dict: Frequencies in a {bitstring: probability, ...} format. + """ + return {bitstring: counts/self.n_shots for bitstring, counts in self.counts.items()} + + def remove_qubit_indices(self, *indices): + """Method to remove qubit indices in the result. The remaining + bitstrings (of new length = length - N_indices) are summed up together. + + Args: + indices (variable number of int): Qubit indices to remove. + """ + new_counts = dict() + for bitstring, counts in self.counts.items(): + new_bitstring = "".join([bitstring[qubit_i] for qubit_i in range(len(bitstring)) if qubit_i not in indices]) + new_counts[new_bitstring] = new_counts.get(new_bitstring, 0) + counts + + self.counts = new_counts + + def post_select(self, expected_outcomes): + """Method to apply post selection on Histogram data, based on a post + selection function in this module. This method changes the Histogram + object inplace. This is explicitly done on desired outcomes for specific + qubit indices. + + Args: + expected_outcomes (dict): Desired outcomes on certain qubit indices + and their expected state. For example, {0: "1", 1: "0"} would + filter results based on the first qubit with the |1> state and + the second qubit with the |0> state measured. + """ + def f_post_select(bitstring): + for qubit_i, expected_bit in expected_outcomes.items(): + if bitstring[qubit_i] != expected_bit: + return False + return True + + new_hist = filter_hist(self, f_post_select) + self.counts = new_hist.counts + self.remove_qubit_indices(*list(expected_outcomes.keys())) + + def resample(self, n_shots): + """Generating new Histogram with n_shots, done with the + get_resampled_frequencies function (see the relevant documentation for + more details). + + Args: + n_shots (int): Self-explanatory. + + Returns: + Histogram: resampled data with n_shots from the distribution. + """ + new_frequencies = get_resampled_frequencies(self.frequencies, n_shots) + return Histogram(new_frequencies, n_shots) + + def get_expectation_value(self, term, coeff=1.): + """Output the expectation value for qubit operator term. The histogram + data is expected to have been acquired in a compatible measurement + basis. + + Args: + term(openfermion-style QubitOperator object): a qubit operator, with + only a single term. + coeff (imaginary): Coefficient to multiply the eigenvalue by. + + Returns: + imaginary: Expectation value for this operator. + """ + return coeff*Simulator.get_expectation_value_from_frequencies_oneterm(term, self.frequencies) + + +def aggregate_histograms(*hists): + """Return a Histogram object formed from data aggregated from all input + Histogram objects. + + Args: + hists (variable number of Histogram objects): the input Histogram objects. + + Returns: + Histogram: The aggregated histogram. + """ + + if not hists: + raise ValueError(f"{aggregate_histograms.__name__} takes at least one Histogram object as an argument") + elif len(hists) == 1: + return hists[0] + else: + # Check that data is on the same number of qubits (same bitstring lengths). + lengths_bitstrings = set([h.n_qubits for h in hists]) + if len(lengths_bitstrings) > 1: + raise ValueError(f"Input histograms have different bitstring lengths ({lengths_bitstrings})") + + # Aggregate histograms. + total_counter = sum([Counter({k: v for k, v in h.counts.items()}) for h in hists], Counter()) + + return Histogram(dict(total_counter)) + + +def filter_hist(hist, function, *args, **kwargs): + """Filter selection function to consider bitstrings in respect to a boolean + predicate on a bitstring. + + Args: + hist (Histogram): Self-explanatory. + function (function): Given a bitstring and some arbitrary arguments, the + predicate should return a boolean value. + + Returns: + Histogram: New Histogram with filtered outcomes based on the predicate + provided. + + """ + new_counts = {bitstring: counts for bitstring, counts in hist.counts.items() if function(bitstring, *args, **kwargs)} + return Histogram(new_counts, msq_first=False) diff --git a/tangelo/toolboxes/post_processing/tests/test_histogram.py b/tangelo/toolboxes/post_processing/tests/test_histogram.py new file mode 100644 index 000000000..c0b2bd26d --- /dev/null +++ b/tangelo/toolboxes/post_processing/tests/test_histogram.py @@ -0,0 +1,200 @@ +# Copyright 2021 Good Chemistry Company. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import unittest + +from tangelo.toolboxes.post_processing import Histogram, aggregate_histograms, filter_hist + + +class HistogramTest(unittest.TestCase): + + def test_init_histogram_from_counts(self): + """Initialize histogram object with a dictionary of counts.""" + h = Histogram({"00": 60, "11": 40}) + + self.assertEqual(h.counts, {"00": 60, "11": 40}) + self.assertEqual(h.n_shots, 100) + self.assertEqual(h.n_qubits, 2) + + ref_freq = {"00": 0.6, "11": 0.4} + for k, v in h.frequencies.items(): + self.assertAlmostEquals(v, ref_freq.get(k)) + + def test_init_histogram_from_frequencies_and_nshots(self): + """Initialize histogram object with a dictonary of probabilities and + a number of shots. + """ + + d = {"00": 0.6, "11": 0.4} + h = Histogram(d, n_shots=100) + + self.assertEqual(h.counts, {"00": 60, "11": 40}) + self.assertEqual(h.n_shots, 100) + self.assertEqual(h.n_qubits, 2) + + ref_freq = {"00": 0.6, "11": 0.4} + for k, v in h.frequencies.items(): + self.assertAlmostEquals(v, ref_freq.get(k)) + + def test_init_histogram_negative_nshots(self): + """Initialize histogram object with a negative n_shots number.""" + + with self.assertRaises(ValueError): + Histogram({"00": 0.6, "11": 0.4}, n_shots=-100) + + def test_histogram_inconsistent_bitstring_length(self): + """Throw an error if bitstrings are not all of consistent length.""" + with self.assertRaises(ValueError): + Histogram({"00": 60, "111": 40}) + + def test_histogram_bitstring_ordering(self): + """Test the behavior with different bit ordering.""" + h = Histogram({"110": 60, "001": 40}, msq_first=True) + self.assertEqual(h.counts, {"011": 60, "100": 40}) + + def test_histogram_bad_frequencies(self): + """Test the behavior of the Histogram initialization with frequencies + not summing up to 1. + """ + with self.assertRaises(ValueError): + Histogram({"00": 0.6, "11": 0.3}, n_shots=100, epsilon=1e-2) + + with self.assertRaises(ValueError): + Histogram({"00": 0.6, "11": 0.5}, n_shots=100, epsilon=1e-2) + + def test_eq_two_histograms(self): + """Test the == operator on two Histograms.""" + d1 = {"00": 60, "11": 40} + d2 = {"00": 50, "11": 50} + + self.assertEqual(Histogram(d1), Histogram(d1)) + self.assertNotEqual(Histogram(d1), Histogram(d2)) + + def test_adding_two_histograms(self): + """Aggregate histograms into a new one.""" + + h1 = Histogram({"00": 60, "11": 40}) + h2 = Histogram({"00": 60, "01": 40}) + + hf = h1 + h2 + self.assertEqual(hf, Histogram({"00": 120, "11": 40, "01": 40})) + + hf += h1 + self.assertEqual(hf, Histogram({"00": 180, "11": 80, "01": 40})) + + def test_removing_one_qubit(self): + """Test removing one bit in an Histogram.""" + + h = Histogram({"00": 40, "01": 30, "10": 20, "11": 10}) + h.remove_qubit_indices(1) + self.assertEqual(h, Histogram({"0": 70, "1": 30})) + + def test_removing_two_qubits(self): + """Test removing two bits in an Histogram.""" + + h = Histogram({"000": 40, "010": 30, "100": 20, "110": 10}) + h.remove_qubit_indices(0, 2) + self.assertEqual(h, Histogram({"0": 60, "1": 40})) + + def test_post_select_method(self): + """Post select in place based on one qubit measurement.""" + + h = Histogram({"00": 40, "01": 30, "10": 20, "11": 10}) + h.post_select(expected_outcomes={0: "0"}) + self.assertEqual(h, Histogram({"0": 40, "1": 30})) + + def test_resample_method(self): + """Test the resampling method. As resampling is probabilistic, we do not + check for the probability outputs. + """ + + h = Histogram({"00": 60, "11": 40}) + + h_10shots = h.resample(10) + + self.assertEqual(h_10shots.n_qubits, 2) + self.assertEqual(h_10shots.n_shots, 10) + + def test_get_expectation_value_method(self): + """Test the output for the get_expectation_value method.""" + + h = Histogram({"00": 50, "11": 50}) + + exp_1z = h.get_expectation_value(((0, "Z"),), 1.) + self.assertAlmostEqual(exp_1z, 0.) + + exp_1zz = h.get_expectation_value(((0, "Z"), (1, "Z")), 1.) + self.assertAlmostEqual(exp_1zz, 1.) + + exp_2zz = h.get_expectation_value(((0, "Z"), (1, "Z")), 2.) + self.assertAlmostEqual(exp_2zz, 2.) + + +class AgregateHistogramsTest(unittest.TestCase): + + def test_aggregate_histograms(self): + """Aggregate histograms into a new one.""" + + h1 = Histogram({"00": 60, "11": 40}) + h2 = Histogram({"00": 60, "01": 40}) + + hf = aggregate_histograms(h1, h2, h1) + self.assertEqual(hf, Histogram({"00": 180, "11": 80, "01": 40})) + + def test_aggregate_histograms_length_inconsistent(self): + """Return an error if histograms have bitstrings of different lengths.""" + with self.assertRaises(ValueError): + h1 = Histogram({"00": 60, "11": 40}) + h2 = Histogram({"1": 100}) + aggregate_histograms(h1, h2) + + def test_aggregate_histograms_empty(self): + """Return an error if no histogram is passed.""" + with self.assertRaises(ValueError): + aggregate_histograms() + + +class FilterTest(unittest.TestCase): + + def test_filter_one_qubit(self): + """Test filter Histogram based on one qubit measurement.""" + + def f(bitstring, qubit_i, outcome): + return bitstring[qubit_i] == outcome + + h = Histogram({"00": 40, "01": 30, "10": 20, "11": 10}) + + h_ps_zero = filter_hist(h, f, 0, "0") + self.assertEqual(h_ps_zero, Histogram({"00": 40, "01": 30})) + + h_ps_one = filter_hist(h, f, 0, "1") + self.assertEqual(h_ps_one, Histogram({"10": 20, "11": 10})) + + def test_filter_many_qubits(self): + """Test filter Histogram based on many qubit measurements.""" + + def f(bitstring, qubit_i, outcome_i, qubit_j, outcome_j): + return (bitstring[qubit_i] == outcome_i) and (bitstring[qubit_j] == outcome_j) + + h = Histogram({"000": 40, "010": 30, "101": 20, "111": 10}) + + h_ps_zero = filter_hist(h, f, 0, "0", 1, "0") + self.assertEqual(h_ps_zero, Histogram({"000": 40})) + + h_ps_one = filter_hist(h, f, 0, "1", 1, "1") + self.assertEqual(h_ps_one, Histogram({"111": 10})) + + +if __name__ == "__main__": + unittest.main()