diff --git a/README.md b/README.md index fbd3e04b5a41aa57142e2449fd59cf8df79dd8bd..4e7c39412c31ef738e64be8cc8c2f4b22290a17c 100644 --- a/README.md +++ b/README.md @@ -40,53 +40,29 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b -p /global/common/software/m4454/env_muraligm source /global/common/software/m4454/env_muraligm/bin/activate - module load cudatoolkit/11.7 + module load cudatoolkit/12.0 ``` - * Make custom conda env - ``` - conda create --name env_custom python=3.9 - conda activate env_custom - ``` - * Install qiskit - ``` - pip install -r requirements.txt --no-cache-dir - ``` - * **NOTE:** Make sure to test the installation with the sample codes provided: - 1. Test qiskit installation: [`test_qiskit_installation.py`](test_qiskit_installation.py) - + * Since the QLSA circuit generator ([HHL](https://github.com/anedumla/quantum_linear_solvers)) uses an older qiskit version requiring qiskit-terra, we need to make two envs: (1) to generate the circuit and (2) to run the circuit using qiskit 1.0 + 1. Install libraries to generate circuit + * Make custom conda env ``` - python test_qiskit_installation.py -backtyp ideal + conda create --name qlsa-circuit python=3.11 + conda activate qlsa-circuit ``` - -
Sample output from the test code: - + * Install qiskit and [linear solver](https://github.com/anedumla/quantum_linear_solvers) package ``` - Backend: QasmSimulator('qasm_simulator') - Job status is JobStatus.DONE - - Total count for 00 and 11 are: {'00': 494, '11': 506} - ┌───┐ ┌─┐ - q_0:┤ H ├──■──┤M├─── - └───┘┌─┴─┐└╥┘┌─┐ - q_1:─────┤ X ├─╫─┤M├ - └───┘ ║ └╥┘ - c:2/═══════════╩══╩═ - 0 1 + pip install -r requirements_circuit.txt --no-cache-dir ``` -
- - * Change `-backtyp` for different backends. - * **NOTE:** To run using IBM Provider, you need to add your IBM Quantum Computing API KEY and instance to the `keys.sh` file and source activate it. - 2. Test [linear solver package](https://github.com/anedumla/quantum_linear_solvers): [`test_linear_solver.py`](test_linear_solver.py) + * Test [linear solver package](https://github.com/anedumla/quantum_linear_solvers): [`test_linear_solver.py`](test_linear_solver.py) ``` python test_linear_solver.py -nq 2 ``` - +
Sample output from the test code: - + ``` Simulator: AerSimulator('aer_simulator') Time elapsed for classical: 0 min 0.00 sec @@ -100,13 +76,13 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an q9_1: ┤1 ├┤5 ├────────┤5 ├ └──────────────┘│ │┌──────┐│ │ q10_0: ───────────────┤0 ├┤3 ├┤0 ├ - │ QPE ││ ││ QPE_dg │ + │ QPE ││ ││ QPE_dg │ q10_1: ───────────────┤1 ├┤2 ├┤1 ├ - │ ││ ││ │ + │ ││ ││ │ q10_2: ───────────────┤2 ├┤1 1/x ├┤2 ├ - │ ││ ││ │ + │ ││ ││ │ q10_3: ───────────────┤3 ├┤0 ├┤3 ├ - └──────┘│ │└─────────┘ + └──────┘│ │└─────────┘ q11: ─────────────── ─────────┤4 ├─────────── └──────┘ tridiagonal state: @@ -116,15 +92,15 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an q82_1:┤1 ├┤5 ├────────┤5 ├ └──────────────┘│ │┌──────┐│ │ q83_0:────────────────┤0 ├┤3 ├┤0 ├ - │ ││ ││ │ + │ ││ ││ │ q83_1:────────────────┤1 QPE ├┤2 ├┤1 QPE_dg ├ - │ ││ ││ │ + │ ││ ││ │ q83_2:────────────────┤2 ├┤1 ├┤2 ├ - │ ││ 1/x ││ │ + │ ││ 1/x ││ │ q83_3:────────────────┤3 ├┤0 ├┤3 ├ - │ ││ ││ │ + │ ││ ││ │ a1: ────────────────┤6 ├┤ ├┤6 ├ - └──────┘│ │└─────────┘ + └──────┘│ │└─────────┘ q84: ────────────────────────┤4 ├─────────── └──────┘ classical Euclidean norm: 1.237833351044751 @@ -143,11 +119,48 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an ===========Data not saved=========== ```
+ + 2. Install libraries to run the circuit + * Make custom conda env + ``` + conda create --name qlsa-solver python=3.11 + conda activate qlsa-solver + ``` + * Install qiskit and other packages + ``` + pip install -r requirements_solver.txt --no-cache-dir + ``` + + * Test qiskit installation: [`test_qiskit_installation.py`](test_qiskit_installation.py) + + ``` + python test_qiskit_installation.py -backtyp ideal + ``` + +
Sample output from the test code: + + ``` + Backend: QasmSimulator('qasm_simulator') + Job status is JobStatus.DONE + + Total count for 00 and 11 are: {'00': 494, '11': 506} + ┌───┐ ┌─┐ + q_0:┤ H ├──■──┤M├─── + └───┘┌─┴─┐└╥┘┌─┐ + q_1:─────┤ X ├─╫─┤M├ + └───┘ ║ └╥┘ + c:2/═══════════╩══╩═ + 0 1 + ``` +
+ + * Change `-backtyp` for different backends. + * **NOTE:** To run using IBM Provider, you need to add your IBM Quantum Computing API KEY and instance to the `keys.sh` file and source activate it. 2. Install GPU version of Aer simulator (skip for Frontier): ``` - pip install qiskit-aer-gpu==0.13.3 --no-cache-dir + pip install qiskit-aer-gpu==0.14.2 --no-cache-dir ```
Notes for Perlmutter: @@ -157,11 +170,15 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an ``` Or ``` - export LD_LIBRARY_PATH=/global/common/software/m4454/env_muraligm/envs/env_custom/lib/python3.9/site-packages/nvidia/cuda_runtime/lib:/global/common/software/m4454/env_muraligm/envs/env_custom/lib/python3.9/site-packages/cuquantum/lib:/global/common/software/m4454/env_muraligm/envs/env_custom/lib/python3.9/site-packages/cutensor/lib:$LD_LIBRARY_PATH + export LD_LIBRARY_PATH=/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/nvidia/cuda_runtime/lib:/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/cuquantum/lib:/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/cutensor/lib:/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/nvidia/nvjitlink/lib:$LD_LIBRARY_PATH ```
- * **NOTE:** Make sure to test the installation with the sample code provided: [`test_gpu.py`](test_gpu.py) + * **NOTE:** Make sure to test the installation: + ``` + python -c "from qiskit_aer import AerSimulator; simulator = AerSimulator(); print(simulator.available_devices())" + ``` + * with the sample code provided: [`test_gpu.py`](test_gpu.py) ``` python test_gpu.py -nq 2 --gpu @@ -234,7 +251,7 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an ``` source /global/common/software/m4454/env_muraligm/bin/activate - module load cudatoolkit/11.7 + module load cudatoolkit/12.0 conda activate env_custom source init_perlmutter.sh ``` @@ -246,9 +263,9 @@ All developments were done on [OLCF Andes](https://docs.olcf.ornl.gov/systems/an ``` * Change `nq` to change size of system of equations -3. Run generalized QLSA script for various use-cases: [`linear_solver.py`](linear_solver.py) +3. Run QLSA circuit generator script for various use-cases: [`circuit_HHL.py`](circuit_HHL.py) ``` - srun -N1 -n1 -c2 python linear_solver.py -case [case-name] -casefile [case-var-file] -s [#-shots] --savedata[optional but recommended] + srun -N1 -n1 -c2 python circuit_HHL.py -case [case-name] -casefile [case-var-file] --savedata ```
Sample output for Hele-Shaw flow, solving for pressure: diff --git a/circuit_HHL.py b/circuit_HHL.py new file mode 100644 index 0000000000000000000000000000000000000000..432a31be4989ea9b4e43501cc09a6e12e888a6c5 --- /dev/null +++ b/circuit_HHL.py @@ -0,0 +1,91 @@ +# Introduction +''' +Script to generate HHL circuit that solves any Ax=b problem. +Function `func_matrix_vector.py` is used to define A and b. +Sample code run script: +python circuit_HHL.py -case sample-tridiag -casefile input_vars.yaml --savedata +''' + +import numpy as np + +# Importing standard Qiskit libraries +from qiskit import QuantumCircuit, transpile +from qiskit import qpy +from qiskit_aer import AerSimulator +from linear_solvers import NumPyLinearSolver, HHL +# library to generate matrix and vector for linear system of equations +import func_matrix_vector as matvec + +import time +import os +import argparse +import pickle + +parser = argparse.ArgumentParser() +parser.add_argument("-case", "--case_name", type=str, default='ideal', required=False, help="Name of the problem case: 'sample-tridiag', 'hele-shaw'") +parser.add_argument("-casefile", "--case_variable_file", type=str, default='ideal', required=False, help="YAML file containing variables for the case: 'input_vars.yaml'") +parser.add_argument("--gpu", default=False, action='store_true', help="Use GPU backend for Aer simulator.") +parser.add_argument("--gpumultiple", default=False, action='store_true', help="Use multiple GPUs for the backend of Aer simulator.") +parser.add_argument("--drawcirc", default=False, action='store_true', help="Draw circuit.") + +parser.add_argument("--savedata", default=False, action='store_true', help="Save data at `models/` with `` based on parameters.") +args = parser.parse_args() + +if __name__ == '__main__': + # Get system matrix and vector + matrix, vector, input_vars = matvec.get_matrix_vector(args) + MATRIX_SIZE = matrix.shape[0] + n_qubits_matrix = int(np.log2(MATRIX_SIZE)) + + # setup quantum backend + backend_type = 'ideal' + backend_method = 'statevector' + print(f'Using \'{backend_type}\' simulator with \'{backend_method}\' backend') + if args.gpu: backend = AerSimulator(method=backend_method, device='GPU') + elif args.gpumultiple: backend = AerSimulator(method=backend_method, device='GPU', blocking_enable=True, blocking_qubits=18) + else: backend = AerSimulator(method=backend_method) + print(f'Backend: {backend}') + + # setup HHL solver + # backend_init = qc_backend('ideal', 'statevector', args) + hhl = HHL(quantum_instance=backend) + + # Solutions + # classical soultion + t = time.time() + classical_solution = NumPyLinearSolver().solve(matrix, vector/np.linalg.norm(vector)) + t_classical = time.time() - t + print(f'Time elapsed for classical: {int(t_classical/60)} min {t_classical%60:.2f} sec', flush=True) + + # generate HHL circuit + print(f'==================Generating HHL circuit================', flush=True) + t = time.time() + circ = hhl.construct_circuit(matrix, vector) + t_circ = time.time() - t + print(f'Time elapsed for generating HHL circuit: {int(t_circ/60)} min {t_circ%60:.2f} sec') + + # Save data + if args.savedata: + circ_transpile = transpile(circ, backend) + # save metadata (DON'T USE Pickle to save the circuit - only works for a given version) + save_data = { 'args' : args, + 'input_vars' : input_vars, + 'matrix' : matrix, + 'vector' : vector, + 't_circ' : t_circ} + filename = input_vars['savefilename'].format(**input_vars) + savefilename = f'{filename}_circ_nqmatrix{n_qubits_matrix}' + file = open(f'{savefilename}.pkl', 'wb') + pickle.dump(save_data, file) + file.close() + # save circuit as QPY file + with open(f'{savefilename}.qpy', 'wb') as fd: + qpy.dump(circ_transpile, fd) + print("===========Circuit saved===========") + + # Plot circuit + if args.drawcirc: + circ.measure_all() + print(f'Circuit:\n{circ.draw()}', flush=True) + + diff --git a/func_matrix_vector.py b/func_matrix_vector.py index 06169d67c30883bbf66cb1921808ccf0fcbb99c0..b90597c7a20bca99cd15df198d4053d377a46276 100644 --- a/func_matrix_vector.py +++ b/func_matrix_vector.py @@ -30,9 +30,9 @@ def sample_tridiag(doc_id, args): input_vars = get_yaml(args.case_variable_file, doc_id) print(f"Case: {input_vars['case_name']}") filename = input_vars['savefilename'].format(**input_vars) - NUM_QUBITS = input_vars['NUM_QUBITS'] + n_qubits_matrix = input_vars['NQ_MATRIX'] # custom systems - MATRIX_SIZE = 2 ** NUM_QUBITS + MATRIX_SIZE = 2 ** n_qubits_matrix # entries of the tridiagonal Toeplitz symmetric matrix a = 1 @@ -41,7 +41,7 @@ def sample_tridiag(doc_id, args): [-1, 0, 1], shape=(MATRIX_SIZE, MATRIX_SIZE)).toarray() vector = np.array([1] + [0]*(MATRIX_SIZE - 1)) - return matrix, vector, filename + return matrix, vector, input_vars def Hele_Shaw(doc_id, args): input_vars = get_yaml(args.case_variable_file, doc_id) @@ -138,7 +138,7 @@ def Hele_Shaw(doc_id, args): if args.savedata == True: MATRIX_SIZE = A_herm.shape[0] - NUM_QUBITS = int(np.log2(MATRIX_SIZE)) + n_qubits_matrix = int(np.log2(MATRIX_SIZE)) save_data = {'P_in' : P_in, 'P_out' : P_out, 'U_top' : U_top, @@ -151,14 +151,14 @@ def Hele_Shaw(doc_id, args): 'ny' : ny, 'A_herm' : A_herm, 'B_herm' : B_herm, - 'NUM_QUBITS' : NUM_QUBITS, + 'n_qubits_matrix' : n_qubits_matrix, 'args' : args} file = open(f'{filename}_metadata.pkl', 'wb') pickle.dump(save_data, file) file.close() print("===========Metadata saved===========") - return A_herm, B_herm, filename + return A_herm, B_herm, input_vars def Cylinder_2D(doc_id, args): input_vars = get_yaml(args.case_variable_file, doc_id) @@ -195,7 +195,7 @@ def Cylinder_2D(doc_id, args): print(f'Reformatted A_herm:\n{A_herm}\nB_herm:\n{B_herm}') print(f'Determinant of resulting matrix: {np.linalg.det(A_herm)}\nCondition # of resulting matrix: {np.linalg.cond(A_herm)}') - return A_herm, B_herm, filename + return A_herm, B_herm, input_vars # Functions def next_power_of_2(x): diff --git a/func_qc.py b/func_qc.py index d20b289c278cd18bdc45ad7d743b4514e6022401..82856fd244cbf394b06c42d9b8f4f982713ec345 100644 --- a/func_qc.py +++ b/func_qc.py @@ -1,23 +1,30 @@ # Introduction ''' -Functions to perform quantum circuit generation for HHL algorith, -transpiling the circuit for tthe specific backend, and +Functions to perform quantum circuit loading, +transpiling the circuit for the specific backend, and running exact and shots-based simulations. +NOTE: The current function qc_circ also computes +the fidelity and solution of a QLSA problem. +Thus, the number of qubits representing the system/matrix +is also needed. If only the output is needed, +any quantum circuit can be loaded and run. ''' import numpy as np -from linear_solvers import NumPyLinearSolver, HHL # Importing standard Qiskit libraries -from qiskit import QuantumCircuit, transpile -from qiskit.execute_function import execute -from qiskit import Aer +from qiskit import QuantumCircuit +from qiskit import qpy from qiskit_aer import AerSimulator +from qiskit_ibm_runtime import SamplerV2 as Sampler +from qiskit_ibm_runtime import RuntimeEncoder +from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager from qiskit.quantum_info import state_fidelity import time import os import argparse import pickle +import json # get backend based on type and method def qc_backend(backend_type, backend_method, args): @@ -28,184 +35,174 @@ def qc_backend(backend_type, backend_method, args): else: backend = AerSimulator(method=backend_method) elif backend_type=='fake': # fake backend - from qiskit.providers import fake_provider - backend = getattr(fake_provider, backend_method)() # FakeNairobi FakePerth FakeMumbai FakeWashington + from qiskit_ibm_runtime.fake_provider import FakeProviderForBackendV2 + backend = FakeProviderForBackendV2().backend(backend_method) elif backend_type=='real-ibm': # real hardware backend - from qiskit_ibm_provider import IBMProvider from qiskit_ibm_runtime import QiskitRuntimeService - # save your IBMProvider accout for future loading + # save your IBM accout for future loading API_KEY = os.getenv('IBMQ_API_KEY') instance = os.getenv('IBMQ_INSTANCE') - IBMProvider.save_account(instance=instance, token=API_KEY, overwrite=True) + # save your QiskitRuntimeService accout for future loading QiskitRuntimeService.save_account( - channel="ibm_quantum", - instance=instance, - token=API_KEY, - overwrite=True + channel="ibm_quantum", + instance=instance, + token=API_KEY, + overwrite=True ) - provider = IBMProvider() # Using IBMProvider to use backend.run() option - backend = provider.get_backend(backend_method) # ibm_nairobi simulator_statevector + service = QiskitRuntimeService() + backend = service.backend(backend_method) else: raise Exception(f'Backend type \'{backend_type}\' not implemented.') return backend # circuit generation, transpile, running -def qc_circ(matrix, vector, hhl, args, backend_method, backend, classical_solution, filename='temp', plot_hist=False): - - MATRIX_SIZE = matrix.shape[0] - NUM_QUBITS = int(np.log2(MATRIX_SIZE)) - print(f'**************************Quantum circuit generation, transpile & running*************************', flush=True) +def qc_circ(n_qubits_matrix, classical_solution, args, input_vars): + ''' + Function to load quantum circuit, transpile, and run. + Input: + n_qubits_matrix num. of quibits representing the system/matrix + classical_solution classical solution of linear system of equations + args input arguments containing details of shots, backend, etc. + input_vars parameters of the system of equations being solved + Output: + job the job handle of Qiskit primitive (only Sampler for now) + ''' + print(f'**************************Quantum circuit loading, transpile & running*************************', flush=True) + # ============================ + # First setup quantum backend + print(f'Using \'{args.backend_type}\' simulator with \'{args.backend_method}\' backend') + backend = qc_backend(args.backend_type, args.backend_method, args) + print(f'Backend: {backend}') + # ============================ - # 1. Generate circuit - savefilename = f'{filename}_circ_nq{NUM_QUBITS}.pkl' - if args.loadcirc == True: - t = time.time() - file = open(savefilename, 'rb') - data = pickle.load(file) - file.close() - circ = data['circ'] - t_load = time.time() - t - print(f'===============Loaded circuit (before transpile) using pickled data==============') - print(f'Time elapsed for loading circuit: {int(t_load/60)} min {t_load%60:.2f} sec', flush=True) - else: - print(f'==================Making a circuit and simulating it================', flush=True) - t = time.time() - circ = hhl.construct_circuit(matrix, vector) - t_circ = time.time() - t - print(f'Time elapsed for generating HHL circuit: {int(t_circ/60)} min {t_circ%60:.2f} sec') - # Save data - if args.savedata == True: - save_data = { 'NUM_QUBITS' : NUM_QUBITS, - 'matrix' : matrix, - 'vector' : vector, - 'circ' : circ, - 't_circ' : t_circ} - file = open(savefilename, 'wb') - pickle.dump(save_data, file) - file.close() - print("===========Circuit saved===========") - # print(circ.qasm()) #filename=f'{savefilename}_qasm') - print(f'Circuit:\n{circ}', flush=True) - if backend_method[0]=='ideal': circ.save_statevector() + # 1. Load generated circuit + filename = input_vars['savefilename'].format(**input_vars) + savefilename = f'{filename}_circ_nqmatrix{n_qubits_matrix}' + t = time.time() + with open(f'{savefilename}.qpy', 'rb') as fd: + circ = qpy.load(fd)[0] + t_load = time.time() - t + print(f'===============Loaded circuit (before transpile) ==============') + print(f'Time elapsed for loading circuit: {int(t_load/60)} min {t_load%60:.2f} sec', flush=True) circ.measure_all() + if args.drawcirc: print(f'Circuit:\n{circ.draw()}', flush=True) + print(f"Circuit details:\n# qubits = {circ.num_qubits}\n# gates = {sum(circ.count_ops().values())}\n# CNOT = {circ.count_ops()['cx']}\nDepth = {circ.depth()}") # ============================ # 2. Transpile circuit for simulator - savefilename = f'{filename}_circ-transpile_nq{NUM_QUBITS}_backend-{backend_method[1]}.pkl' - if args.loadcirctranspile == True: + savefilename = f'{filename}_circ-transpile_nqmatrix{n_qubits_matrix}_backend-{args.backend_method}' + if args.loadcirctranspile: t = time.time() - file = open(savefilename, 'rb') - data = pickle.load(file) - file.close() - circ = data['circ'] + with open(f'{savefilename}.qpy', 'rb') as fd: + isa_circ = qpy.load(fd)[0] t_load = time.time() - t - print(f'===============Loaded transpiled circuit using pickled data==============') + print(f'===============Loaded transpiled circuit ==============') print(f'Time elapsed for loading circuit: {int(t_load/60)} min {t_load%60:.2f} sec', flush=True) else: t = time.time() - circ = transpile(circ, backend) + # Convert to ISA circuits: Circuits must obey the ISA of the backend. + pm = generate_preset_pass_manager(backend=backend, optimization_level=1) + isa_circ = pm.run(circ) t_transpile = time.time() - t print(f'Time elapsed for transpiling the circuit: {int(t_transpile/60)} min {t_transpile%60:.2f} sec') # Save data - if args.savedata == True: - save_data = { 'NUM_QUBITS' : NUM_QUBITS, - 'matrix' : matrix, - 'vector' : vector, - 'circ' : circ, - 't_circ' : t_circ, + if args.savedata: + save_data = { 'args' : args, + 'input_vars' : input_vars, 't_transpile' : t_transpile} - file = open(savefilename, 'wb') + file = open(f'{savefilename}.pkl', 'wb') pickle.dump(save_data, file) file.close() + # save transpiled circuit + with open(f'{savefilename}.qpy', 'wb') as fd: + qpy.dump(isa_circ, fd) print("===========Transpiled Circuit saved===========", flush=True) # ============================ # 3. Run and get counts shots = args.SHOTS + # setup Sampler + sampler = Sampler(backend) + t = time.time() - result = backend.run(circ, shots=shots, memory=True).result() - # result = execute(circ, backend, shots=shots, memory=True).result() - # extracting the counts for the given number of counts based on the probabilities obtained from the true simulation + # Run the job + job = sampler.run([isa_circ]) + # Grab results from the job + result = job.result() t_run = time.time() - t print(f'Time elapsed for running the circuit: {int(t_run/60)} min {t_run%60:.2f} sec', flush=True) - counts = result.get_counts(circ) + # Returns counts + counts = result[0].data.meas.get_counts() print(f'counts:\n{counts}') - # Returning measurement outcomes for each shot - memory = result.get_memory(circ) - # print(f'memory: {memory}') - # Saving the final statevector if using ideal (qiskit) backend - if backend_method[0]=='ideal': - # statevector = result.get_statevector(circ) # remove circ to get results of shots-based simulation - statevector = result.get_statevector() + if args.backend_type=='ideal': + isa_circ.remove_final_measurements() # no measurements allowed + from qiskit.quantum_info import Statevector + statevector = Statevector(isa_circ) statevector = np.asarray(statevector) istart = int(len(statevector)/2) - exact_vector = statevector[istart:istart+MATRIX_SIZE].real + exact_vector = statevector[istart:istart+(int(2**n_qubits_matrix))].real # get counts based probabilistic/statistical state vector - counts_ancilla, counts_total, probs_vector, counts_vector = get_ancillaqubit(counts, NUM_QUBITS) + counts_ancilla, counts_total, probs_vector, counts_vector = get_ancillaqubit(counts, n_qubits_matrix) print(f'All counts of ancila (only the first 2**nq represent solution vector):\n{counts_ancilla}') print("Counts vector should approach exact vector in infinite limit") print(f'counts_vector:\n{counts_vector}') - if backend_method[0]=='ideal': print(f'exact_vector/norm:\n{exact_vector/np.linalg.norm(exact_vector)}') + if args.backend_type=='ideal': print(f'exact_vector/norm:\n{exact_vector/np.linalg.norm(exact_vector)}') # print solutions - print(f'\ntrue solution:\n{classical_solution.state}') + print(f'\ntrue solution:\n{classical_solution}') # normalize counts vector with true solution norm - counts_solution_vector = classical_solution.euclidean_norm * counts_vector / np.linalg.norm(counts_vector) + counts_solution_vector = np.linalg.norm(classical_solution) * counts_vector / np.linalg.norm(counts_vector) print(f'\ncounts solution vector:\n{counts_solution_vector}') - print(f'diff with true solution (%):\n{np.abs(classical_solution.state-counts_solution_vector)*100/(classical_solution.state+1e-15)}') - print(f'Fidelity: {fidelity(counts_solution_vector, classical_solution.state)}') - if backend_method[0]=='ideal': - exact_solution_vector = classical_solution.euclidean_norm * exact_vector / np.linalg.norm(exact_vector) + print(f'diff with true solution (%):\n{np.abs(classical_solution-counts_solution_vector)*100/(classical_solution+1e-15)}') + print(f'Fidelity: {fidelity(counts_solution_vector, classical_solution)}') + if args.backend_type=='ideal': + exact_solution_vector = np.linalg.norm(classical_solution) * exact_vector / np.linalg.norm(exact_vector) print(f'\nexact solution vector:\n{exact_solution_vector}') - print(f'diff with true solution (%):\n{np.abs(classical_solution.state-exact_solution_vector)*100/(classical_solution.state+1e-15)}') - print(f'Fidelity: {fidelity(exact_solution_vector, classical_solution.state)}') + print(f'diff with true solution (%):\n{np.abs(classical_solution-exact_solution_vector)*100/(classical_solution+1e-15)}') + print(f'Fidelity: {fidelity(exact_solution_vector, classical_solution)}') # plot histogram - if plot_hist: - from qiskit.tools.visualization import plot_histogram + if args.plothist: + from qiskit.visualization import plot_histogram import matplotlib.pyplot as plt - plot_histogram(counts, figsize=(7, 7), color='tab:green', title=f'{backend_method[0]}:{backend_method[1]}') # dodgerblue tab:green + plot_histogram(counts, figsize=(7, 7), color='tab:green', title=f'{args.backend_type}:{args.backend_method}') # dodgerblue tab:green plt.savefig('Figs/temp_hist.png') # Save full data - savefilename = f'{filename}_circ-fullresults_nq{NUM_QUBITS}_backend-{backend_method[1]}_shots{shots}.pkl' - if args.savedata == True: - save_data = { 'NUM_QUBITS' : NUM_QUBITS, - 'matrix' : matrix, - 'vector' : vector, - 'circ' : circ, - 'shots' : shots, - 'result' : result, + savefilename = f'{filename}_circ-fullresults_nqmatrix{n_qubits_matrix}_backend-{args.backend_method}_shots{shots}' + if args.savedata: + save_data = { 'args' : args, + 'input_vars' : input_vars, 'counts' : counts, - 'memory' : memory, 'exact_vector' : exact_vector, 'counts_ancilla' : counts_ancilla, 'counts_vector' : counts_vector, 'counts_solution_vector' : counts_solution_vector, - 'exact_solution_vector' : exact_solution_vector, 'classical_solution' : classical_solution, - 't_circ' : t_circ, - 't_transpile' : t_transpile, 't_run' : t_run} - file = open(savefilename, 'wb') - pickle.dump(save_data, file) - file.close() + with open(f"{savefilename}.pkl", "wb") as file: + pickle.dump(save_data, file) + # save results + with open(f"{savefilename}_result.json", "w") as file: + json.dump(result, file, cls=RuntimeEncoder) print("===========Full data saved===========") + return job + # function to measure the qubits def get_ancillaqubit(counts, nq): ''' NOTE: only count measurements when ancilla qubit (leftmost) is 1 - input: + Input: counts counts from the simulator nq number of qubits used to represent the system or solution vector - output: + Output: counts_ancill acounts of the measurements where ancilla qubit = 1 other metricis for combination of nq qubits = 1 ''' @@ -240,6 +237,14 @@ def get_ancillaqubit(counts, nq): # function to compute fidelity of the solution def fidelity(qfunc, true): + ''' + Function to compute fidelity of solution state. + Input: + qfunc quantum solution state + true classiccal/true solution state + Output: + fidelity state fidelity + ''' solution_qfun_normed = qfunc / np.linalg.norm(qfunc) solution_true_normed = true / np.linalg.norm(true) fidelity = state_fidelity(solution_qfun_normed, solution_true_normed) diff --git a/init_perlmutter.sh b/init_perlmutter.sh index 010f45f23364775121940fa9b809738f8189b75a..28d4d5c700874dab0c4072bf0e320ddc41fe3de1 100644 --- a/init_perlmutter.sh +++ b/init_perlmutter.sh @@ -1,5 +1,4 @@ #!/bin/bash -export LD_LIBRARY_PATH=/global/common/software/m4454/env_muraligm/envs/env_custom/lib/python3.9/site-packages/nvidia/cuda_runtime/lib:/global/common/software/m4454/env_muraligm/envs/env_custom/lib/python3.9/site-packages/cuquantum/lib:/global/common/software/m4454/env_muraligm/envs/env_custom/lib/python3.9/site-packages/cutensor/lib:$LD_LIBRARY_PATH - +export LD_LIBRARY_PATH=/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/nvidia/cuda_runtime/lib:/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/cuquantum/lib:/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/cutensor/lib:/global/common/software/m4454/env_muraligm/envs/qlsa-solver/lib/python3.11/site-packages/nvidia/nvjitlink/lib:$LD_LIBRARY_PATH diff --git a/input_vars.yaml b/input_vars.yaml index 8e7061e4fe503769703eda8d2026dc5bb731ed94..27063252e99cddf203d7548659cee804bdf3991f 100644 --- a/input_vars.yaml +++ b/input_vars.yaml @@ -1,7 +1,7 @@ --- # Case 0: Sample system - tridiagonal matrix case_name: Sample tridiagonal -NUM_QUBITS: 2 # Numer of qubits to determine size of linear system of quations being solved. A matrix size = 2^NUM_QUBITS. +NQ_MATRIX: 2 # Numer of qubits to determine size of linear system of quations being solved. A matrix size = 2^NQ_MATRIX. savedir: models savefilename: "{savedir}/sample_HHL" diff --git a/requirements.txt b/requirements_circuit.txt similarity index 100% rename from requirements.txt rename to requirements_circuit.txt diff --git a/requirements_solver.txt b/requirements_solver.txt new file mode 100644 index 0000000000000000000000000000000000000000..8f2f17e10cc3f875a74c3146fd565a4dfb7901d7 --- /dev/null +++ b/requirements_solver.txt @@ -0,0 +1,10 @@ +PyYAML==6.0.1 +qiskit==1.0.2 +qiskit-aer==0.14.1 +qiskit-ibm-runtime==0.23.0 +qiskit-qasm3-import==0.4.2 +jupyterlab==4.2.0 +matplotlib==3.8.4 + + + diff --git a/linear_solver.py b/solver.py similarity index 54% rename from linear_solver.py rename to solver.py index 0db9a9abb8ed5325144fa3844026aa420dbd6c37..84473e68ab9c919fff6d52f69184a924de44ecec 100644 --- a/linear_solver.py +++ b/solver.py @@ -4,18 +4,13 @@ Script to perform quantum linear solver for any Ax=b problem. The HHL algorithm Function `func_matrix_vector.py` is used to define A and b. Function `func_qc.py` is used to generate, transpile, and run the quantum circuit. Sample code run script: -python linear_solver.py -case sample-tridiag -casefile input_vars.yaml -s 1000 +python solver.py -case sample-tridiag -casefile input_vars.yaml -s 1000 ''' import numpy as np # Importing standard Qiskit libraries -from qiskit import QuantumCircuit, transpile -from qiskit.execute_function import execute -from qiskit import Aer -from qiskit_aer import AerSimulator -from linear_solvers import NumPyLinearSolver, HHL -# library to generate HHL circuit for given matrix and vector, transpile circuit with given backend, and run shot-based simulations +# library to load circuit for given matrix and vector, transpile circuit with given backend, and run shot-based simulations from func_qc import qc_backend, qc_circ import func_matrix_vector as matvec @@ -32,36 +27,28 @@ parser.add_argument("-s", "--SHOTS", type=int, default=1000, required=True, help parser.add_argument("--gpu", default=False, action='store_true', help="Use GPU backend for Aer simulator.") parser.add_argument("--gpumultiple", default=False, action='store_true', help="Use multiple GPUs for the backend of Aer simulator.") parser.add_argument("-backtyp", "--backend_type", type=str, default='ideal', required=False, help="Type of the backend: 'ideal', 'fake' 'real-ibm'") -parser.add_argument("-backmet", "--backend_method", type=str, default='statevector', required=False, help="Method/name of the backend. E.g. 'statevector', 'simulator_statevector' 'FakeNairobi' 'ibm_nairobi'") +parser.add_argument("-backmet", "--backend_method", type=str, default='statevector', required=False, help="Method/name of the backend. E.g. 'statevector' 'fake_sherbrooke' 'ibm_sherbrooke'") +parser.add_argument("--drawcirc", default=False, action='store_true', help="Draw circuit.") +parser.add_argument("--plothist", default=False, action='store_true', help="Draw circuit.") parser.add_argument("--savedata", default=False, action='store_true', help="Save data at `models/` with `` based on parameters.") -parser.add_argument("--loadcirc", default=False, action='store_true', help="Load circuit at `models/` with `` based on parameters.") parser.add_argument("--loadcirctranspile", default=False, action='store_true', help="Load transpiled circuit at `models/` with `` based on parameters.") args = parser.parse_args() -# Get system matrix and vector -matrix, vector, filename = matvec.get_matrix_vector(args) - -# setup HHL solver -backend_init = qc_backend('ideal', 'statevector', args) -hhl = HHL(quantum_instance=backend_init) - -# setup quantum backend -backend_type = args.backend_type -backend_method = args.backend_method -print(f'Using \'{backend_type}\' simulator with \'{backend_method}\' backend') -backend = qc_backend(backend_type, backend_method, args) -print(f'Backend: {backend}') - -# Solutions -# classical soultion -t = time.time() -classical_solution = NumPyLinearSolver().solve(matrix, vector/np.linalg.norm(vector)) -t_classical = time.time() - t -print(f'Time elapsed for classical: {int(t_classical/60)} min {t_classical%60:.2f} sec', flush=True) - -# quantum circuit solution -qc_circ(matrix, vector, hhl, args, [backend_type, backend_method], backend, classical_solution, filename=filename) +if __name__ == '__main__': + # Get system matrix and vector + matrix, vector, input_vars = matvec.get_matrix_vector(args) + n_qubits_matrix = int(np.log2(matrix.shape[0])) + + # Solutions + # classical soultion + t = time.time() + classical_solution = np.linalg.solve(matrix, vector/np.linalg.norm(vector)) + t_classical = time.time() - t + print(f'Time elapsed for classical: {int(t_classical/60)} min {t_classical%60:.2f} sec', flush=True) + + # quantum solution + qc_circ(n_qubits_matrix, classical_solution, args, input_vars) diff --git a/test_gpu.py b/test_gpu.py index e12b4b96a7930399d22a8b3a353086be15e28bb9..56985da3b5969ee2149ef0fe9178bcb48246841b 100644 --- a/test_gpu.py +++ b/test_gpu.py @@ -3,10 +3,10 @@ Sample code to test the GPU (and multi-GPU) capability of aer_simulator Reference: https://qiskit.org/ecosystem/aer/howtos/running_gpu.html ''' -from qiskit import QuantumCircuit, transpile +from qiskit import QuantumCircuit from qiskit_aer import AerSimulator from qiskit.circuit.library import QuantumVolume -from qiskit.execute_function import execute +from qiskit import transpile import argparse import time @@ -16,22 +16,23 @@ parser.add_argument("--gpu", default=False, action='store_true', help="Use GPU b parser.add_argument("--gpumultiple", default=False, action='store_true', help="Use multiple GPUs for the backend of Aer simulator.") args = parser.parse_args() -if args.gpu: sim = AerSimulator(method='statevector', device='GPU') -elif args.gpumultiple: sim = AerSimulator(method='statevector', device='GPU', blocking_enable=True, blocking_qubits=18) -else: sim = AerSimulator(method='statevector') -print(f'Simulator: {sim}') +# Select backend +if args.gpu: backend = AerSimulator(method='statevector', device='GPU') +elif args.gpumultiple: backend = AerSimulator(method='statevector', device='GPU', blocking_enable=True, blocking_qubits=18) +else: backend = AerSimulator(method='statevector') +print(f'Simulator: {backend}') +# Generate circuit and transpile qubits = args.NUM_QUBITS t = time.time() -circ = transpile(QuantumVolume(qubits, seed = 0)) -circ.measure_all() +qc = QuantumVolume(qubits, seed = 0) +qc.measure_all() +qc = transpile(qc, backend) elpsdt1 = time.time() - t - +# Run circuit t = time.time() -# if args.gpu: result = execute(circ, sim, shots=100, blocking_enable=True, blocking_qubits=23).result() -# else: result = execute(circ, sim, shots=100).result() -result = execute(circ, sim, shots=100).result() +result = backend.run(qc).result() elpsdt2 = time.time() - t print(f'N qubits: {qubits}; GPU: {args.gpu}; multiple-GPU: {args.gpumultiple};\nTime elapsed 1: {int(elpsdt1/60)} min {elpsdt1%60:.2f} sec\nTime elapsed 2: {int(elpsdt2/60)} min {elpsdt2%60:.2f} sec') diff --git a/test_linear_solver.py b/test_linear_solver.py index 071af91ead518fa82a53c2d9c6fb1107b9b0b78d..192a5c2edb38dc3e86c70f7e690f9c979e4f90f2 100644 --- a/test_linear_solver.py +++ b/test_linear_solver.py @@ -7,7 +7,7 @@ import numpy as np from scipy.sparse import diags # Importing Qiskit libraries -from qiskit import QuantumCircuit, transpile +from qiskit import QuantumCircuit from qiskit_aer import AerSimulator from linear_solvers import NumPyLinearSolver, HHL from linear_solvers.matrices.tridiagonal_toeplitz import TridiagonalToeplitz @@ -18,7 +18,7 @@ import argparse import pickle parser = argparse.ArgumentParser() -parser.add_argument("-nq", "--NUM_QUBITS", type=int, default=2, required=True, help="Numer of qubits to determine size of linear system of quations (A*x=b) being solved. Size of A matrix = 2^NUM_QUBITS.") +parser.add_argument("-nq", "--NQ_MATRIX", type=int, default=2, required=True, help="Numer of qubits to determine size of linear system of quations (A*x=b) being solved. Size of A matrix = 2^NQ_MATRIX.") parser.add_argument("--savedata", default=False, action='store_true', help="Save data at `models/` with `` based on parameters.") parser.add_argument("--gpu", default=False, action='store_true', help="Use GPU backend for Aer simulator.") args = parser.parse_args() @@ -32,8 +32,8 @@ args = parser.parse_args() # vector = np.array([1, 0]) # tridi_matrix = TridiagonalToeplitz(1, 1, -1 / 3) # custom systems -NUM_QUBITS = args.NUM_QUBITS -MATRIX_SIZE = 2 ** NUM_QUBITS +n_qubits_matrix = args.NQ_MATRIX +MATRIX_SIZE = 2 ** n_qubits_matrix # entries of the tridiagonal Toeplitz symmetric matrix a = 1 b = -1/3 @@ -42,14 +42,11 @@ matrix = diags([b, a, b], shape=(MATRIX_SIZE, MATRIX_SIZE)).toarray() vector = np.array([1] + [0]*(MATRIX_SIZE - 1)) # we also generate an optimized matrix construction - tridiagonal toeplitz -tridi_matrix = TridiagonalToeplitz(NUM_QUBITS, a, b) +tridi_matrix = TridiagonalToeplitz(n_qubits_matrix, a, b) # ============ # Select backend: Using different simulators (default in `linear_solvers` is statevector simulation) -from qiskit import Aer -# https://qiskit.org/documentation/tutorials/simulators/1_aer_provider.html -# run `Aer.backends()` to see simulators -backend = Aer.get_backend('aer_simulator_statevector') +backend = AerSimulator(method='statevector') if args.gpu: backend.set_options(device='GPU') backend.set_options(precision='single') @@ -146,7 +143,7 @@ print(f'diff (%): {np.abs(classical_solution.state-solvec_tridi)*100/classical_s savedata = args.savedata if savedata == True: savedir = f'models' - savefilename = f'{savedir}/sample_HHL_numq{NUM_QUBITS}_fulldata' + savefilename = f'{savedir}/sample_HHL_numq{n_qubits_matrix}_fulldata' if not os.path.exists(savedir): os.mkdir(savedir) n=2 @@ -155,7 +152,7 @@ if savedata == True: n+=1 savefilename += '.pkl' - save_data = { 'NUM_QUBITS' : NUM_QUBITS, + save_data = { 'n_qubits_matrix' : n_qubits_matrix, 'a' : a, 'b' : b, 'matrix' : matrix, diff --git a/test_qiskit_installation.py b/test_qiskit_installation.py index 47417a65f8c26bff652a9b1e84915a4eeff8011b..033d1caac055bc1367cf5c2db0ecc51829d41874 100644 --- a/test_qiskit_installation.py +++ b/test_qiskit_installation.py @@ -1,8 +1,9 @@ import numpy as np import time import os -from qiskit import QuantumCircuit, transpile -from qiskit.providers.jobstatus import JobStatus +from qiskit import QuantumCircuit +from qiskit_ibm_runtime import SamplerV2 as Sampler +from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager import argparse parser = argparse.ArgumentParser() @@ -12,33 +13,31 @@ args = parser.parse_args() backend_type = args.backend_type # Choose the simulator or backend to run the quantum circuit if backend_type=='ideal': - # Using ideal simulator, Aer's qasm_simulator (works even without IBMQ account, don't have to wait in a queue) - from qiskit.providers.aer import QasmSimulator - backend = QasmSimulator() + # Using ideal simulator, AerSimulator (works even without IBMQ account, don't have to wait in a queue) + from qiskit_aer import AerSimulator + backend = AerSimulator() elif backend_type=='fake': - # Using qiskit's fake provider (works even without IBMQ account, don't have to wait in a queue) - from qiskit.providers import fake_provider - backend = getattr(fake_provider, 'FakeNairobi')() # FakePerth FakeMumbai FakeWashington + # Using qiskit's fake provider (works even without IBMQ account, don't have to wait in a queue) + from qiskit_ibm_runtime.fake_provider import FakeProviderForBackendV2 # or use FakeProvider for V1 fake backends - https://docs.quantum.ibm.com/api/qiskit-ibm-runtime/fake_provider + backend = FakeProviderForBackendV2().backend("fake_sherbrooke") elif backend_type=='real-ibm': - # Using the latest qiskit_ibm_provider - #### IF YOU HAVE AN IBMQ ACCOUNT (using an actual backend) ##### - - from qiskit_ibm_provider import IBMProvider - from qiskit_ibm_runtime import QiskitRuntimeService - # save your IBMProvider accout for future loading - API_KEY = os.getenv('IBMQ_API_KEY') - instance = os.getenv('IBMQ_INSTANCE') - IBMProvider.save_account(instance=instance, token=API_KEY, overwrite=True) - - # save your QiskitRuntimeService accout for future loading - QiskitRuntimeService.save_account( + # Using the latest qiskit_ibm_runtime + #### IF YOU HAVE AN IBMQ ACCOUNT (using an actual backend) ##### + + from qiskit_ibm_runtime import QiskitRuntimeService + # save your IBM accout for future loading + API_KEY = os.getenv('IBMQ_API_KEY') + instance = os.getenv('IBMQ_INSTANCE') + + # save your QiskitRuntimeService accout for future loading + QiskitRuntimeService.save_account( channel="ibm_quantum", instance=instance, token=API_KEY, overwrite=True - ) - provider = IBMProvider() - backend = provider.get_backend("simulator_statevector") # ibm_nairobi simulator_statevector + ) + service = QiskitRuntimeService() + backend = service.backend("ibm_sherbrooke") else: raise Exception(f'Backend type \'{backend_type}\' not implemented.') @@ -46,7 +45,7 @@ print(f'Backend: {backend}') ###################################### # Create a Quantum Circuit acting on the q register -circuit = QuantumCircuit(2, 2) +circuit = QuantumCircuit(2) # Add a H gate on qubit 0 circuit.h(0) @@ -55,28 +54,21 @@ circuit.h(0) circuit.cx(0, 1) # Map the quantum measurement to the classical bits -circuit.measure([0,1], [0,1]) +circuit.measure_all() -# compile the circuit down to low-level QASM instructions -# supported by the backend (not needed for simple circuits) -compiled_circuit = transpile(circuit, backend) - -# Execute the circuit on the qasm simulator -job = backend.run(compiled_circuit, shots=1000) - -# Make a "waiting in queue" message -while job.status() is not JobStatus.DONE: - print("Job status is", job.status() ) - time.sleep(30) - -print("Job status is", job.status() ) +# Circuits must obey the ISA of the backend. +# Convert to ISA circuits +pm = generate_preset_pass_manager(backend=backend, optimization_level=1) +isa_circuit = pm.run(circuit) +# setup Sampler +sampler = Sampler(backend) +# Run the job +job = sampler.run([isa_circuit]) # Grab results from the job result = job.result() - # Returns counts -counts = result.get_counts(compiled_circuit) -print("\nTotal count for 00 and 11 are:",counts) +print(f"\nTotal count for 00 and 11 are: {result[0].data.meas.get_counts()}") # Draw the circuit print(circuit.draw())