ZNE is incompatible with AdaptVQE()

[MarcoBarroca]

Describe the bug
I have been trying to use AdaptVQE() with ZNE by setting the resilience_level but have been running into the following error:

2023-03-30T15:28:36.439337986Z Traceback (most recent call last):
2023-03-30T15:28:36.439366661Z   File "/provider/programruntime/program_starter_wrapper.py", line 91, in execute
2023-03-30T15:28:36.439416418Z     final_result = self.main(backend, self.messenger, **self.user_params)
2023-03-30T15:28:36.439427471Z   File "/code/program.py", line 1507, in main
2023-03-30T15:28:36.439438964Z     result = job.result()
2023-03-30T15:28:36.439450175Z   File "/opt/app-root/lib64/python3.9/site-packages/zne/meta/job.py", line 62, in result
2023-03-30T15:28:36.439460209Z     result = self.zne_strategy.mitigate_noisy_result(result)
2023-03-30T15:28:36.439469620Z   File "/opt/app-root/lib64/python3.9/site-packages/zne/zne_strategy.py", line 257, in mitigate_noisy_result
2023-03-30T15:28:36.439481389Z     val, meta = self.extrapolator.extrapolate_zero(data)
2023-03-30T15:28:36.439493970Z   File "/opt/app-root/lib64/python3.9/site-packages/zne/extrapolation/extrapolator.py", line 78, in extrapolate_zero
2023-03-30T15:28:36.439532518Z     return self.infer(0, data)
2023-03-30T15:28:36.439543991Z   File "/opt/app-root/lib64/python3.9/site-packages/zne/extrapolation/extrapolator.py", line 70, in infer
2023-03-30T15:28:36.439554151Z     model, metadata = self.fit_regression_model(data)
2023-03-30T15:28:36.439562787Z   File "/opt/app-root/lib64/python3.9/site-packages/zne/extrapolation/extrapolator.py", line 56, in fit_regression_model
2023-03-30T15:28:36.439572911Z     data = self._validate_data(data)
2023-03-30T15:28:36.439587111Z   File "/opt/app-root/lib64/python3.9/site-packages/zne/extrapolation/extrapolator.py", line 102, in _validate_data
2023-03-30T15:28:36.439597188Z     raise TypeError(
2023-03-30T15:28:36.439637867Z TypeError: Invalid data provided, expeceted Sequence[Sequence[float, float, float]].

This seems related to the way that the results from the job are formatted as the error happens during the extrapolation procedure while it tries to read the results.

Steps to reproduce
I'll attach a python script in the next comment to reproduce the issue. Note that it requires PySCF and qiskit-nature.

The retry_primitves.py is so that it can implement a max_retries option on the estimator and was taken from #598. The problem persists even if you use the standard Runtime Estimator().

Expected behavior
Ideally, it would work like it does when you use standard VQE. You can modify the attached script to use Standard VQE and it runs without issues.

Suggested solutions
I'm unsure if the modification would need to be in AdaptVQE() or in the ZNE files. I've taken a look at the files in the zne-prototype repo hoping they would be similar but it didn't help to come up with a solution.

Additional Information

• qiskit-ibm-runtime version: 0.9.1
• Python version: 3.9
• Operating system: Mac OS 12.6 arm64

[MarcoBarroca]

Python script to reproduce the issue

from qiskit.algorithms import VQE, NumPyMinimumEigensolver, NumPyEigensolver #Algorithms

#Qiskit odds and ends
from qiskit.circuit.library import EfficientSU2, EvolvedOperatorAnsatz
from qiskit.algorithms.optimizers import COBYLA, SPSA, SLSQP, L_BFGS_B
from qiskit.opflow import Z2Symmetries, X, Y, Z, I, PauliSumOp, Gradient, NaturalGradient
from qiskit import IBMQ, BasicAer, Aer, transpile
from qiskit.utils import QuantumInstance, algorithm_globals
from qiskit.utils.mitigation import CompleteMeasFitter #Measurement error mitigatioin
from qiskit.tools.visualization import circuit_drawer
from qiskit.providers.aer import AerSimulator
from qiskit.providers.aer.noise import NoiseModel
from qiskit.algorithms.minimum_eigensolvers import VQE, AdaptVQE, MinimumEigensolverResult
from qiskit.primitives import Estimator
from qiskit_aer.primitives import Estimator as AerEstimator
from qiskit.quantum_info import SparsePauliOp

#qiskit_nature
from qiskit_nature.second_q.drivers import PySCFDriver, MethodType
from qiskit_nature.second_q.formats.molecule_info import MoleculeInfo
from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.circuit.library import UCCSD, PUCCD, SUCCD, HartreeFock, CHC, VSCF
from qiskit_nature.second_q.operators.fermionic_op import FermionicOp
from qiskit_nature.second_q.transformers import ActiveSpaceTransformer , FreezeCoreTransformer
from qiskit_nature.second_q.problems import ElectronicStructureProblem, EigenstateResult
from qiskit_nature.second_q.mappers import QubitConverter, ParityMapper, BravyiKitaevMapper, JordanWignerMapper
from qiskit_nature.second_q.algorithms.ground_state_solvers.minimum_eigensolver_factories.vqe_ucc_factory import VQEUCCFactory
from qiskit_nature.second_q.algorithms.ground_state_solvers.minimum_eigensolver_factories.numpy_minimum_eigensolver_factory import NumPyMinimumEigensolverFactory
from qiskit_nature.second_q.algorithms.ground_state_solvers import GroundStateEigensolver
from qiskit_nature.second_q.algorithms.excited_states_solvers.eigensolver_factories.numpy_eigensolver_factory import NumPyEigensolverFactory
from qiskit_nature.second_q.algorithms.excited_states_solvers import QEOM, ExcitedStatesEigensolver

#Runtime
from qiskit_ibm_runtime import (QiskitRuntimeService, Session,
                                Estimator as RuntimeEstimator)
from qiskit_ibm_runtime.options import Options, ResilienceOptions, SimulatorOptions, TranspilationOptions, ExecutionOptions 
from retry_primitives import RetryEstimator

#PySCF
from functools import reduce
import scipy.linalg
from pyscf import scf
from pyscf import gto, dft
from pyscf import mcscf, fci
from functools import reduce
from pyscf.mcscf import avas, dmet_cas

#Python odds and ends
import pylab
import numpy as np
import os
from IPython.display import display, clear_output
import mapomatic as mm

from datetime import datetime

## Python program to store list to file using pickle module
import pickle

# write list to binary file
def write_list(a_list,filename):
    # store list in binary file so 'wb' mode
    with open(filename, 'wb') as fp:
        pickle.dump(a_list, fp)
        print('Done writing list into a binary file')
        
def write_dict(a_dict,filename):
    # store list in binary file so 'wb' mode
    with open(filename, 'wb') as fp:
        pickle.dump(a_dict, fp,protocol=pickle.HIGHEST_PROTOCOL)
        print('Done writing dict into a binary file')

# Read list to memory
def read(filename):
    # for reading also binary mode is important
    with open(filename, 'rb') as fp:
        n_list = pickle.load(fp)
        return n_list

IBMQ.load_account()
service = QiskitRuntimeService(channel='ibm_quantum')
seed=42

def callback(eval_count, param, val,meta):  
    # Overwrites the same line when printing
    counts.append(eval_count)
    interim_info['counts'].append(eval_count)
    values.append(val)
    interim_info['values'].append(val)
    params.append(param)
    interim_info['params'].append(param)
    meta_dicts.append(meta)
    mean=np.mean(values)
    std=np.std(values)

    write_dict(interim_info,f'MG+H2O_mont_4qubit_zne_interim_info_{int(distances[0]*10)}')
    display("Evaluation: {}, Energy: {}, Mean: {}, STD: {}, Metadata: {}".format(eval_count, val,mean,std, meta))
    clear_output(wait=True)

def adapt_solver(distances,mapper,optimizer,freeze_core):
    
    dists=[]
    results=[]
    problems=[]
    hf_energies=[]
    initial_points=[]
    ops=[]
    ansatze=[]
    for dist in distances:
        #Driver
        driver,og_problem=make_driver(dist)
        #Qubit_Op
        qubit_op, problem, converter,hf_energy = make_qubit_op(dist,og_problem,mapper,freeze_core)
        ops.append(qubit_op)
        problems.append(problem)
        #Initial State
        init_state = HartreeFock(problem.num_spatial_orbitals, problem.num_particles, converter)
        
        #ansatz
        ansatz = UCCSD(num_spatial_orbitals=problem.num_spatial_orbitals,num_particles=problem.num_particles,qubit_converter=converter,initial_state=init_state)

        operator_pool = []
        for op in ansatz.operators:
            sop = op.primitive
            for pauli, coeff in zip(sop.paulis, sop.coeffs):
                if sum(pauli.x & pauli.z) % 2 == 0:
                    continue
            operator_pool.append(PauliSumOp(coeff * SparsePauliOp(pauli)))

        ansatz = EvolvedOperatorAnsatz(
        operators=operator_pool,
        initial_state=init_state,
        )
        ansatze.append(ansatz)
        
        
        # Set initial parameters of the ansatz
        if len(initial_points)!=0 and ops[-1].num_qubits==ops[-2].num_qubits:
            initial_point=initial_points[-1]
        elif len(initial_points)!=0:
            old_ans=ansatz[-1]
            if ansatz.num_parameters>old_ans.num_parameters:
                initial_point=np.append(initial_points[-1],np.zeros(ansatz.num_parameters - old_ans.num_parameters))
            else:
                to_remove=old_ans.num_parameters-ansatz.num_parameters
                initial_point=np.delete(initial_points[-1],np.arange(-1,-to_remove-1,-1))
        else:
            #initial_point = np.pi/4 * np.random.rand(ansatz.num_parameters)
            initial_point= np.zeros(ansatz.num_parameters)
        
        estimator = Estimator([ansatz], [qubit_op])
    
        counts = []
        values = []
        deviation = []
    
        custom_vqe =VQE(estimator=estimator, ansatz=ansatz, optimizer=optimizer,initial_point=initial_point, callback=callback)
        
        adapt_vqe=AdaptVQE(custom_vqe)
        adapt_vqe.supports_aux_operators = lambda: True
        adapt_vqe.threshold=1e-3
        adapt_vqe.max_iterations=4
        adapt_vqe.delta=1e-4
        
        solver = GroundStateEigensolver(converter, adapt_vqe)
        result = solver.solve(problem)
        
        initial_points.append(result.raw_result.optimal_point)
        results.append(result)
        dists.append(dist)
        hf_energies.append(hf_energy)
    
    return results, problems ,distances, hf_energies


def real_solver(distances, mapper, optimizer,freeze_core,est_options, device):
    
    dists=[]
    results=[]
    problems=[]
    hf_energies=[]
    initial_points=[]
    ops=[]
    ansatze=[]
                             
    for dist in distances:
        #Driver
        driver,og_problem=make_driver(dist)
        #Qubit_Op
        qubit_op, problem, converter,hf_energy = make_qubit_op(dist,og_problem,mapper,freeze_core)
        ops.append(qubit_op)
        problems.append(problem)
        #Initial State
        init_state = HartreeFock(problem.num_spatial_orbitals, problem.num_particles, converter)
        
        #ansatz
        #ansatz = UCCSD(num_spatial_orbitals,num_particles,converter)
        
        ansatz = UCCSD(num_spatial_orbitals=problem.num_spatial_orbitals,num_particles=problem.num_particles,qubit_converter=converter,initial_state=init_state)

        operator_pool = []
        for op in ansatz.operators:
            sop = op.primitive
            for pauli, coeff in zip(sop.paulis, sop.coeffs):
                if sum(pauli.x & pauli.z) % 2 == 0:
                    continue
            operator_pool.append(PauliSumOp(coeff * SparsePauliOp(pauli)))

        ansatz = EvolvedOperatorAnsatz(
        operators=operator_pool,
        initial_state=init_state,
        )
        ansatze.append(ansatz)
        #ansatz = EfficientSU2(num_qubits=qubit_op.num_qubits,su2_gates='ry', entanglement='linear', reps=3, initial_state=init_state)
        
        #ansatz_opt = transpile(ansatz, backend=provider.get_backend(device),optimization_level=3,routing_method='sabre')
        #small_qc = mm.deflate_circuit(ansatz_opt)
        #layouts = mm.matching_layouts(small_qc, provider.get_backend(device))
        #scores = mm.evaluate_layouts(small_qc, layouts, provider.get_backend(device))
        #ansatz = transpile(small_qc, backend=provider.get_backend(device),initial_layout=scores[0][0],optimization_level=3,routing_method='sabre')
        
        # Set initial parameters of the ansatz
        if len(initial_points)!=0 and ops[-1].num_qubits==ops[-2].num_qubits:
            initial_point=initial_points[-1]
        elif len(initial_points)!=0:
            old_ans=ansatze[-1]
            if ansatz.num_parameters>old_ans.num_parameters:
                initial_point=np.append(initial_points[-1],np.zeros(ansatz.num_parameters - old_ans.num_parameters))
            else:
                to_remove=old_ans.num_parameters-ansatz.num_parameters
                initial_point=np.delete(initial_points[-1],np.arange(-1,-to_remove-1,-1))
        else:
            #initial_point = np.pi/4 * np.random.rand(ansatz.num_parameters)
            initial_point= np.zeros(ansatz.num_parameters)
    
        counts= []
        values = []
        deviation = []
        params=[]
    
        with Session(service=service, backend=device, max_time=288000) as session:
            # Prepare primitive
            rt_estimator = RetryEstimator(session=session,options=est_options,max_retries=6,timeout=18000)
            # Set up algorithm
            custom_vqe = VQE(rt_estimator, ansatz, optimizer,initial_point=initial_point, callback=callback)
            adapt_vqe=AdaptVQE(custom_vqe)
            adapt_vqe.supports_aux_operators = lambda: True
            adapt_vqe.threshold=1e-3
            adapt_vqe.max_iterations=4
            adapt_vqe.delta=1e-4
            # Run algorithm
            solver = GroundStateEigensolver(converter, adapt_vqe)
            result = solver.solve(problem)
            #result = custom_vqe.compute_minimum_eigenvalue(qubit_op,initial_point)
            
            
        initial_points.append(result.raw_result.optimal_point)
        results.append(result)
        dists.append(dist)
        hf_energies.append(hf_energy)
        
        problems.append(problem)
       
    return results,problems, distances

def make_driver(d):
    
    
    molecule = MoleculeInfo(
             # coordinates in Angstrom
                     symbols=['O','H','H','Mg'],
                     coords=[
                            # (d+0.504284,0.0,0.758602),
                            # (d,0.0,0.0),
                            # (d+2*0.504284,0.0,0.0),
                            # (0.0, 0.0, 0.0),
                            (0.0,0.0,0.0),
                            (-0.504284,0.0,-0.758602),
                            (0.504284,0.0,-0.758602),
                            (0.0, 0.0, d),
                            ],
                     multiplicity=1,  # = 2*spin + 1
                     charge=2,
                     units=DistanceUnit.ANGSTROM
                    )
    
    #Set driver
    #driver = PySCFDriver.from_molecule(molecule, basis="sto3g", method=MethodType.ROHF)
    #driver.xc_functional='pbe,pbe'
    driver = PySCFDriver.from_molecule(molecule, basis="6-31g*", method=MethodType.ROKS)
    driver.xc_functional='b3lyp'
    driver.conv_tol = 1e-6

    #Get properties
    problem = driver.run()
    

    return driver, problem

def make_qubit_op(d,og_problem, mapper,freeze_core):
    mol = gto.Mole()
    mol.atom = [
        # ['O',(d+0.504284,0.0,0.758602)],
        # ['H',(d,0.0,0.0),],
        # ['H',(d+2*0.504284,0.0,0.0)],
        # ['Mg',(0.0, 0.0, 0.0)]
        ['O',(0.0,0.0,0)],
        ['H',(-0.504284,0.0,-0.758602),],
        ['H',(0.504284,0.0,-0.758602)],
        ['Mg',(0.0, 0.0, d)]
        ]
    mol.charge=2
    mol.basis = '6-31g*'
    mol.spin = 0
    mol.build()
    
    #mf= scf.ROHF(mol).x2c()
    mf = dft.ROKS(mol).density_fit(auxbasis='def2-universal-jfit')
    mf.xc ='pbe,pbe'
    mf.max_cycle = 50
    mf.conv_tol = 1e-6
    
    first_run=mf.kernel()
    a = mf.stability()[0]
    if(mf.converged):
        energy=first_run
    else:
        mf.max_cycle = 80
        mf.conv_tol = 1e-6
        mf = scf.newton(mf)
        scnd_run=mf.kernel(dm0 = mf.make_rdm1(a,mf.mo_occ)) # using rdm1 constructed from stability analysis
      #mf.kernel(mf.make_rdm1()) #using the rdm from the non-converged calculation
        if(mf.converged):
            energy=scnd_run
        else:
            mf.conv_tol = 1e-6
            mf.max_cycle = 80
            mf = scf.newton(mf) #Second order solver
            energy=mf.kernel(dm0 = mf.make_rdm1())


    ao_labels = ['Mg 2p', 'O 2p']
    avas_obj = avas.AVAS(mf, ao_labels)
    avas_obj.kernel()
    weights=np.append(avas_obj.occ_weights,avas_obj.vir_weights)
    weights=(weights>0.2)*weights
    orbs=np.nonzero(weights)
    orbs=np.nonzero(weights)
    
    transformer = ActiveSpaceTransformer(
            num_electrons=(3,3), #Electrons in active space
            num_spatial_orbitals=4, #Orbitals in active space
            #active_orbitals=orbs[0].tolist().append(orbs[0][-1]+1)
        )
    fz_transformer=FreezeCoreTransformer(freeze_core=freeze_core)
    
    #Define the problem

    problem=transformer.transform(og_problem)
    if freeze_core==True:
        problem=fz_transformer.transform(problem)
        converter = QubitConverter(mapper)
    else:
        converter = QubitConverter(mapper,two_qubit_reduction=True, z2symmetry_reduction='auto')

    hamiltonian=problem.hamiltonian
    second_q_op = hamiltonian.second_q_op()
    
    num_spatial_orbitals = problem.num_spatial_orbitals
    num_particles = problem.num_particles
    
    qubit_op = converter.convert(second_q_op,num_particles=num_particles,sector_locator=problem.symmetry_sector_locator)
    
        

    
    return qubit_op, problem,  converter, energy
    
    
#Runtime Estimator options
ts_opt=TranspilationOptions(
    skip_transpilation=False
)

res_opt=ResilienceOptions(
    #noise_factors=tuple(range(1, 6, 2)),
    noise_amplifier='LocalFoldingAmplifier'
    #extrapolator='LinearExtrapolator'
)

ex_opt=ExecutionOptions(
    shots=1024
)
est_options=Options(
    resilience_level=2,
    optimization_level=3,
    execution=ex_opt,
    #resilience=res_opt,
    #transpilation=ts_opt
)


distances = [0.3,0.6,0.9,1.2,1.5,1.8,2.1,2.4,2.7,3.0,3.3]
optimizer = SPSA(maxiter=100)
mapper=ParityMapper()

algorithm_globals.random_seed = seed

counts = []
values = []
params = []
meta_dicts=[]
interim_info={'counts':[],
              'values':[],
              'params':[]
             }


vqe_results,vqe_problems,dists=real_solver(distances=[distances[0]],
                                            mapper=mapper,
                                            optimizer=optimizer,
                                            freeze_core=False,
                                            est_options=est_options,
                                            device='ibmq_montreal'
                                          )

write_list(vqe_results,f'MG+H2O_mont_4qubit_zne_vqe_results_{int(distances[0]*10)}')
write_list(vqe_problems,f'MG+H2O_mont_4qubit_zne_vqe_problems_{int(distances[0]*10)}')

retry_primitives.py:

import signal, time

from qiskit_ibm_runtime import Sampler, Estimator, Session
from qiskit.providers import JobStatus

def timeout_handler(signum, frame):
    raise Exception('Iteration timed out')

class RetryPrimitiveMixin:
    """RetryPrimitive class.
    
    This class inherits from Qiskit IBM Runtime's Primitives and overwrites its run method such that it retries calling it
    a maximum of 'max_retries' consecutive times, if it encounters one of the following randomly occuring errors:
    
    * A Primitive error (in this case "Job.ERROR" is printed, and the job is cancelled automatically)
    * A timeout error where the job either remains running or completes but does not return anything, for a time larger 
      than 'timeout' (in this case the job is cancelled by the patch and "Job.CANCELLED" is printed)
    * A creation error, where the job fails to be created because connection is lost between the runtime server and the
      quantum computer (in this case "Failed to create job." is printed). If this error occurs, the patch connects the user
      to a new Session (to be handled with care! also, this will unfortunately put the next job in the queue). 
    """
    
    def __init__(self, *args, max_retries: int = 5, timeout: int = 3600, **kwargs) -> None:
        super().__init__(*args, **kwargs)
        self.max_retries = max_retries
        self.timeout = timeout
        self.backend = super().session._backend
        signal.signal(signal.SIGALRM, timeout_handler)
        
    def run(self, *args, **kwargs):
        result = None
        for i in range(self.max_retries):
            try:
                job = super().run(*args, **kwargs)
                while job.status() in [JobStatus.INITIALIZING, JobStatus.QUEUED, JobStatus.VALIDATING]:
                    time.sleep(5) # Check every 5 seconds whether job status has changed
                signal.alarm(self.timeout) # Once job starts running, set timeout to 1 hour by default
                result = job.result()
                if result is not None:
                    signal.alarm(0) # Reset timer
                    return job
            except Exception as e:
                signal.alarm(0) # Reset timer
                print("\nSomething went wrong...")
                print(f"\n\nERROR MESSAGE:\n{e}\n\n")
                if 'job' in locals(): # Sometimes job fails to create
                    print(f"Job ID: {job.job_id}. Job status: {job.status()}.")
                    if job.status() not in [JobStatus.DONE, JobStatus.ERROR, JobStatus.CANCELLED]:
                        job.cancel()
                else:
                    print("Failed to create job.")
                    try:
                        super().session.close()
                        print("Current session was closed.")
                    except:
                        print("Current session could not be closed. Will leave it to close automatically.")
                    print(f"Creating new session...\n")
                    self._session = Session(backend=self.backend)
                print(f"Starting trial number {i+2}...\n")
                signal.alarm(0) # Reset timer
        if result is None:
            raise RuntimeError(f"Program failed! Maximum number of retries ({self.max_retries}) exceeded")
            
class RetrySampler(RetryPrimitiveMixin, Sampler):
    pass

class RetryEstimator(RetryPrimitiveMixin, Estimator):
    pass

[jyu00]

transferred to qiskit-algorithms: qiskit-community/qiskit-algorithms#191 since it looks like AdaptVQE is not passing the correct extrapolators