MarcoBarroca
commented
1 year ago 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
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:
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