get noise working

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 = <unsigned char>(probs[jj] < rand_vals[nprobs*ii+jj])
+            temp = <unsigned char>(probs[jj] > rand_vals[nprobs*ii+jj])
             if temp:
                 mem[ii,mem_slots[jj]] = temp

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 = <unsigned char>(probs[jj] < rand_vals[nprobs*ii+jj])
+            temp = <unsigned char>(probs[jj] > rand_vals[nprobs*ii+jj])
             if temp:
                 mem[ii,memory_slots[jj]] = temp

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
 

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.

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

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

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, '<string>', '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,41 +129,47 @@ 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.")
 
                 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
                 _p = np.cumsum(n_dp / np.sum(n_dp))
                 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
+        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

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,18 +103,38 @@ 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,
                                    task_args=(self.op_system.global_data,
                                    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

qiskit/providers/aer/openpulse/solver/unitary.py

@@ -64,6 +64,11 @@ def unitary_evolution(exp, global_data, ode_options):
     code = compile(_inst, '<string>', '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:]: