From 208ddd735ea8f7ddf790f0fb995e9cbefe80b6b3 Mon Sep 17 00:00:00 2001 From: Paul Nation Date: Tue, 11 Jun 2019 13:14:44 -0400 Subject: [PATCH] get noise working --- .../aer/openpulse/cython/measure.pyx | 2 +- .../providers/aer/openpulse/cython/memory.pyx | 2 +- qiskit/providers/aer/openpulse/qobj/digest.py | 6 +- .../providers/aer/openpulse/qobj/operators.py | 2 +- .../providers/aer/openpulse/solver/codegen.py | 11 +- .../aer/openpulse/solver/data_config.py | 13 ++- .../aer/openpulse/solver/monte_carlo.py | 105 ++++++++++-------- .../providers/aer/openpulse/solver/opsolve.py | 31 +++++- .../providers/aer/openpulse/solver/unitary.py | 5 + 9 files changed, 115 insertions(+), 62 deletions(-) diff --git a/qiskit/providers/aer/openpulse/cython/measure.pyx b/qiskit/providers/aer/openpulse/cython/measure.pyx index 3191d87e01..0c5fda6be1 100644 --- a/qiskit/providers/aer/openpulse/cython/measure.pyx +++ b/qiskit/providers/aer/openpulse/cython/measure.pyx @@ -59,6 +59,6 @@ def write_shots_memory(unsigned char[:, ::1] mem, cdef unsigned char temp for ii in range(nrows): for jj in range(nprobs): - temp = (probs[jj] < rand_vals[nprobs*ii+jj]) + temp = (probs[jj] > rand_vals[nprobs*ii+jj]) if temp: mem[ii,mem_slots[jj]] = temp \ No newline at end of file diff --git a/qiskit/providers/aer/openpulse/cython/memory.pyx b/qiskit/providers/aer/openpulse/cython/memory.pyx index 6cbfbd9d38..fe762ca543 100644 --- a/qiskit/providers/aer/openpulse/cython/memory.pyx +++ b/qiskit/providers/aer/openpulse/cython/memory.pyx @@ -39,6 +39,6 @@ def write_memory(unsigned char[:, ::1] mem, cdef unsigned char temp for ii in range(nrows): for jj in range(nprobs): - temp = (probs[jj] < rand_vals[nprobs*ii+jj]) + temp = (probs[jj] > rand_vals[nprobs*ii+jj]) if temp: mem[ii,memory_slots[jj]] = temp \ No newline at end of file diff --git a/qiskit/providers/aer/openpulse/qobj/digest.py b/qiskit/providers/aer/openpulse/qobj/digest.py index fe59a4791c..1ccace4545 100644 --- a/qiskit/providers/aer/openpulse/qobj/digest.py +++ b/qiskit/providers/aer/openpulse/qobj/digest.py @@ -44,7 +44,7 @@ def digest_pulse_obj(qobj): # Look for config keys out.global_data['shots'] = 1024 if 'shots' in config_keys: - out.global_data['shots'] = qobj['config']['shots'] + out.global_data['shots'] = int(qobj['config']['shots']) out.global_data['seed'] = None if 'seed' in config_keys: @@ -123,7 +123,9 @@ def digest_pulse_obj(qobj): noise.parse() out.noise = noise.compiled - out.can_sample = False + if any(out.noise): + out.can_sample = False + out.global_data['c_num'] = len(out.noise) else: out.noise = None diff --git a/qiskit/providers/aer/openpulse/qobj/operators.py b/qiskit/providers/aer/openpulse/qobj/operators.py index 2edb9985e9..b180128c45 100644 --- a/qiskit/providers/aer/openpulse/qobj/operators.py +++ b/qiskit/providers/aer/openpulse/qobj/operators.py @@ -65,7 +65,7 @@ def gen_oper(opname, index, h_osc, h_qub, states=None): return op.tensor(opers) -def qubit_occ_oper(target_qubit, h_osc, h_qub, level=1): +def qubit_occ_oper(target_qubit, h_osc, h_qub, level=0): """Builds the occupation number operator for a target qubit in a qubit oscillator system, where the oscillator are the first subsystems, and the qubit last. diff --git a/qiskit/providers/aer/openpulse/solver/codegen.py b/qiskit/providers/aer/openpulse/solver/codegen.py index f3e73600a9..4c472b29a1 100644 --- a/qiskit/providers/aer/openpulse/solver/codegen.py +++ b/qiskit/providers/aer/openpulse/solver/codegen.py @@ -40,7 +40,7 @@ def __init__(self, op_system): self.op_system = op_system self.dt = op_system.dt - self.num_ham_terms = len(self.op_system.system) + self.num_ham_terms = self.op_system.global_data['num_h_terms'] # Code generator properties self._file = None @@ -157,6 +157,15 @@ def func_vars(self): spmv_str = "spmvpy(&data{i}[0], &idx{i}[0], &ptr{i}[0], "+ \ "&vec[0], 1.0, &out[0], num_rows)" func_vars.append(spmv_str.format(i=idx)) + + # There is a noise term + if len(self.op_system.system) < self.num_ham_terms: + spmv_str = "spmvpy(&data{i}[0], &idx{i}[0], &ptr{i}[0], "+ \ + "&vec[0], 1.0, &out[0], num_rows)" + func_vars.append("") + func_vars.append("# Noise term") + func_vars.append(spmv_str.format(i=idx+1)) + return func_vars diff --git a/qiskit/providers/aer/openpulse/solver/data_config.py b/qiskit/providers/aer/openpulse/solver/data_config.py index 87fea5aec8..d526fa8fb7 100644 --- a/qiskit/providers/aer/openpulse/solver/data_config.py +++ b/qiskit/providers/aer/openpulse/solver/data_config.py @@ -31,15 +31,17 @@ def op_data_config(op_system): op_system.global_data['c_num'] = 0 if op_system.noise: op_system.global_data['c_num'] = len(op_system.noise) + op_system.global_data['num_h_terms'] += 1 op_system.global_data['c_ops_data'] = [] op_system.global_data['c_ops_ind'] = [] op_system.global_data['c_ops_ptr'] = [] op_system.global_data['n_ops_data'] = [] op_system.global_data['n_ops_ind'] = [] - op_system.global_data['n_ops_ind'] = [] + op_system.global_data['n_ops_ptr'] = [] # if there are any collapse operators + H_noise = 0 for kk in range(op_system.global_data['c_num']): c_op = op_system.noise[kk] n_op = c_op.dag() * c_op @@ -50,10 +52,13 @@ def op_data_config(op_system): # norm ops op_system.global_data['n_ops_data'].append(n_op.data.data) op_system.global_data['n_ops_ind'].append(n_op.data.indices) - op_system.global_data['n_ops_ind'].append(n_op.data.indptr) + op_system.global_data['n_ops_ptr'].append(n_op.data.indptr) # Norm ops added to time-independent part of # Hamiltonian to decrease norm - H[0] -= 0.5j * n_op + H_noise -= 0.5j * n_op + + if H_noise: + H = H + [H_noise] # construct data sets op_system.global_data['h_ops_data'] = [-1.0j* hpart.data.data for hpart in H] @@ -63,7 +68,7 @@ def op_data_config(op_system): # setup ode args string ode_var_str = "" # Hamiltonian data - for kk in range(num_h_terms): + for kk in range(op_system.global_data['num_h_terms']): h_str = "global_data['h_ops_data'][%s], " % kk h_str += "global_data['h_ops_ind'][%s], " % kk h_str += "global_data['h_ops_ptr'][%s], " % kk diff --git a/qiskit/providers/aer/openpulse/solver/monte_carlo.py b/qiskit/providers/aer/openpulse/solver/monte_carlo.py index bc1d5cf8a7..cf1c56f70f 100644 --- a/qiskit/providers/aer/openpulse/solver/monte_carlo.py +++ b/qiskit/providers/aer/openpulse/solver/monte_carlo.py @@ -18,46 +18,51 @@ # Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson. # All rights reserved. +from math import log import numpy as np from scipy.integrate import ode from scipy.linalg.blas import get_blas_funcs from qutip.cy.spmatfuncs import cy_expect_psi_csr, spmv, spmv_csr from openpulse.solver.zvode import qiskit_zvode +from openpulse.cython.memory import write_memory +from openpulse.cython.measure import occ_probabilities, write_shots_memory dznrm2 = get_blas_funcs("znrm2", dtype=np.float64) -def monte_carlo(pid, ophandler, ode_options): +def monte_carlo(seed, exp, global_data, ode_options): """ Monte Carlo algorithm returning state-vector or expectation values at times tlist for a single trajectory. """ - global _cy_rhs_func - - tlist = ophandler._data['tlist'] - memory = [ - [ - 0 for _l in range(ophandler.qobj_config['memory_slot_size']) - ] for _r in range(ophandler.qobj_config['memory_slots']) - ] - register = bytearray(ophandler.backend_config['n_registers']) - - opt = ophandler._options + cy_rhs_func = global_data['rhs_func'] + rng = np.random.RandomState(seed) + tlist = exp['tlist'] + snapshots = [] + # Init memory + memory = np.zeros((1, global_data['memory_slots']), dtype=np.uint8) + # Init register + register = np.zeros(global_data['n_registers'], dtype=np.uint8) + + # Get number of acquire, snapshots, and conditionals + num_acq = len(exp['acquire']) + acq_idx = 0 + num_snap = len(exp['snapshot']) + snap_idx = 0 + num_cond = len(exp['cond']) + cond_idx = 0 collapse_times = [] collapse_operators = [] - # SEED AND RNG AND GENERATE - prng = RandomState(ophandler._options.seeds[pid]) # first rand is collapse norm, second is which operator - rand_vals = prng.rand(2) + rand_vals = rng.rand(2) - ODE = ode(_cy_rhs_func) + ODE = ode(cy_rhs_func) - _inst = 'ODE.set_f_params(%s)' % ophandler._data['string'] + _inst = 'ODE.set_f_params(%s)' % global_data['string'] code = compile(_inst, '', 'exec') exec(code) - psi = ophandler._data['psi0'] # initialize ODE solver for RHS ODE._integrator = qiskit_zvode(method=ode_options.method, @@ -69,27 +74,27 @@ def monte_carlo(pid, ophandler, ode_options): min_step=ode_options.min_step, max_step=ode_options.max_step ) - + # Forces complex ODE solving if not len(ODE._y): ODE.t = 0.0 ODE._y = np.array([0.0], complex) ODE._integrator.reset(len(ODE._y), ODE.jac is not None) + + ODE.set_initial_value(global_data['initial_state'], 0) # make array for collapse operator inds - cinds = np.arange(ophandler._data['c_num']) - n_dp = np.zeros(ophandler._data['c_num'], dtype=float) + cinds = np.arange(global_data['c_num']) + n_dp = np.zeros(global_data['c_num'], dtype=float) - kk_prev = 0 # RUN ODE UNTIL EACH TIME IN TLIST - for kk in tlist: - ODE.set_initial_value(psi, kk_prev) + for stop_time in tlist: # ODE WHILE LOOP FOR INTEGRATE UP TO TIME TLIST[k] - while ODE.t < kk: + while ODE.t < stop_time: t_prev = ODE.t y_prev = ODE.y norm2_prev = dznrm2(ODE._y) ** 2 - # integrate up to kk, one step at a time. - ODE.integrate(kk, step=1) + # integrate up to stop_time, one step at a time. + ODE.integrate(stop_time, step=1) if not ODE.successful(): raise Exception("ZVODE step failed!") norm2_psi = dznrm2(ODE._y) ** 2 @@ -99,11 +104,11 @@ def monte_carlo(pid, ophandler, ode_options): # ------------------------------------------------ ii = 0 t_final = ODE.t - while ii < ophandler._options.norm_steps: + while ii < ode_options.norm_steps: ii += 1 t_guess = t_prev + \ - math.log(norm2_prev / rand_vals[0]) / \ - math.log(norm2_prev / norm2_psi) * (t_final - t_prev) + log(norm2_prev / rand_vals[0]) / \ + log(norm2_prev / norm2_psi) * (t_final - t_prev) ODE._y = y_prev ODE.t = t_prev ODE._integrator.call_args[3] = 1 @@ -113,7 +118,7 @@ def monte_carlo(pid, ophandler, ode_options): "ZVODE failed after adjusting step size!") norm2_guess = dznrm2(ODE._y)**2 if (abs(rand_vals[0] - norm2_guess) < - ophandler._options.norm_tol * rand_vals[0]): + ode_options.norm_tol * rand_vals[0]): break elif (norm2_guess < rand_vals[0]): # t_guess is still > t_jump @@ -124,7 +129,7 @@ def monte_carlo(pid, ophandler, ode_options): t_prev = t_guess y_prev = ODE.y norm2_prev = norm2_guess - if ii > ophandler._options.norm_steps: + if ii > ode_options.norm_steps: raise Exception("Norm tolerance not reached. " + "Increase accuracy of ODE solver or " + "Options.norm_steps.") @@ -132,9 +137,9 @@ def monte_carlo(pid, ophandler, ode_options): collapse_times.append(ODE.t) # all constant collapse operators. for i in range(n_dp.shape[0]): - n_dp[i] = cy_expect_psi_csr(ophandler._data['n_ops_data'][i], - ophandler._data['n_ops_ind'][i], - ophandler._data['n_ops_ptr'][i], + n_dp[i] = cy_expect_psi_csr(global_data['n_ops_data'][i], + global_data['n_ops_ind'][i], + global_data['n_ops_ptr'][i], ODE._y, 1) # determine which operator does collapse and store it @@ -142,23 +147,29 @@ def monte_carlo(pid, ophandler, ode_options): j = cinds[_p >= rand_vals[1]][0] collapse_operators.append(j) - state = spmv_csr(ophandler._data['c_ops_data'][j], - ophandler._data['c_ops_ind'][j], - ophandler._data['c_ops_ptr'][j], + state = spmv_csr(global_data['c_ops_data'][j], + global_data['c_ops_ind'][j], + global_data['c_ops_ptr'][j], ODE._y) state /= dznrm2(state) ODE._y = state ODE._integrator.call_args[3] = 1 - rand_vals = prng.rand(2) + rand_vals = rng.rand(2) # after while loop (Do measurement or conditional) - # ---------------- + # ------------------------------------------------ out_psi = ODE._y / dznrm2(ODE._y) - - # measurement - psi = _proj_measurement(pid, ophandler, kk, out_psi, memory, register) - - kk_prev = kk - - return psi, memory \ No newline at end of file + for aind in range(acq_idx, num_acq): + if exp['acquire'][aind][0] == stop_time: + current_acq = exp['acquire'][aind] + qubits = current_acq[1] + memory_slots = current_acq[2] + probs = occ_probabilities(qubits, out_psi, global_data['measurement_ops']) + rand_vals = rng.rand(memory_slots.shape[0]) + write_shots_memory(memory, memory_slots, probs, rand_vals) + int_mem = memory.dot(np.power(2.0, + np.arange(memory.shape[1]-1,-1,-1))).astype(int) + acq_idx += 1 + + return int_mem \ No newline at end of file diff --git a/qiskit/providers/aer/openpulse/solver/opsolve.py b/qiskit/providers/aer/openpulse/solver/opsolve.py index a21c707bb0..28c4de3a55 100644 --- a/qiskit/providers/aer/openpulse/solver/opsolve.py +++ b/qiskit/providers/aer/openpulse/solver/opsolve.py @@ -32,6 +32,7 @@ from openpulse.solver.data_config import op_data_config import openpulse.solver.settings as settings from openpulse.solver.unitary import unitary_evolution +from openpulse.solver.monte_carlo import monte_carlo from qiskit.tools.parallel import parallel_map dznrm2 = get_blas_funcs("znrm2", dtype=np.float64) @@ -53,12 +54,12 @@ def opsolve(op_system): if not op_system.ode_options.num_cpus: op_system.ode_options.num_cpus = qutip.settings.num_cpus + # build Hamiltonian data structures + op_data_config(op_system) # compile Cython RHS _op_generate_rhs(op_system) # Load cython function _op_func_load(op_system) - # build Hamiltonian data structures - op_data_config(op_system) # load monte carlo class mc = OP_mcwf(op_system) # Run the simulation @@ -102,10 +103,10 @@ def __init__(self, op_system): def run(self): + map_kwargs = {'num_processes': self.op_system.ode_options.num_cpus} + # If no collapse terms, and only measurements at end # can do a single shot. - map_kwargs = {'num_processes': self.op_system.ode_options.num_cpus} - if self.op_system.can_sample: results = parallel_map(unitary_evolution, self.op_system.experiments, @@ -113,7 +114,27 @@ def run(self): self.op_system.ode_options ), **map_kwargs ) - + + # need to simulate each trajectory, so shots*len(experiments) times + # Do a for-loop over experiments, and do shots in parallel_map + else: + results = [] + for exp in self.op_system.experiments: + rng = np.random.RandomState(exp['seed']) + seeds = rng.randint(np.iinfo(np.int32).max-1, + size=self.op_system.global_data['shots']) + exp_res = parallel_map(monte_carlo, + seeds, + task_args=(exp, self.op_system.global_data, + self.op_system.ode_options + ), **map_kwargs + ) + unique = np.unique(exp_res, return_counts=True) + hex_dict = {} + for kk in range(unique[0].shape[0]): + key = hex(unique[0][kk]) + hex_dict[key] = unique[1][kk] + results.append(hex_dict) _cython_build_cleanup(self.op_system.global_data['rhs_file_name']) return results diff --git a/qiskit/providers/aer/openpulse/solver/unitary.py b/qiskit/providers/aer/openpulse/solver/unitary.py index 68e2d851f8..fc84784c07 100644 --- a/qiskit/providers/aer/openpulse/solver/unitary.py +++ b/qiskit/providers/aer/openpulse/solver/unitary.py @@ -64,6 +64,11 @@ def unitary_evolution(exp, global_data, ode_options): code = compile(_inst, '', 'exec') exec(code) + if not len(ODE._y): + ODE.t = 0.0 + ODE._y = np.array([0.0], complex) + ODE._integrator.reset(len(ODE._y), ODE.jac is not None) + # Since all experiments are defined to start at zero time. ODE.set_initial_value(global_data['initial_state'], 0) for kk in tlist[1:]: