Import backend; load circuit; transpile using pm; run using Sampler

# 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

import qiskit.qasm3

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):

        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",

          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):

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}')

    MATRIX_SIZE = matrix.shape[0]

    NUM_QUBITS = int(np.log2(MATRIX_SIZE))

    print(f'**************************Quantum circuit generation, transpile & running*************************', flush=True)

    # ============================

    # 1. Generate circuit

    savefilename = f'{filename}_circ_nq{NUM_QUBITS}.pkl'

    if args.loadcirc == True:

    # 1. Load generated circuit

    filename = input_vars['savefilename'].format(**input_vars)

    savefilename = f'{filename}_circ_nqmatrix{n_qubits_matrix}'

    t = time.time()

        file = open(savefilename, 'rb')

        data = pickle.load(file)

        file.close()

        circ = data['circ']

    circ = qiskit.qasm3.load(f'{savefilename}.qasm')

    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()

    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']

        isa_circ = qiskit.qasm3.load(f'{savefilename}.qasm')

        t_load = time.time() - t

        print(f'===============Loaded transpiled circuit using pickled data==============')

        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}.qasm', 'w') as ofile:

                qiskit.qasm3.dump(isa_circ, ofile)

            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')

        with open(f"{savefilename}.pkl", "wb") as file:

            pickle.dump(save_data, file)

        file.close()

        # 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

    '''

# 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)