diff --git a/pygsti/baseobjs/basis.py b/pygsti/baseobjs/basis.py
index 5141de9d3..08af5324a 100644
--- a/pygsti/baseobjs/basis.py
+++ b/pygsti/baseobjs/basis.py
@@ -333,7 +333,7 @@ def elsize(self):
         int
         """
         if self.elshape is None: return 0
-        return int(_np.product(self.elshape))
+        return int(_np.prod(self.elshape))
 
     @property
     def first_element_is_identity(self):
@@ -957,7 +957,7 @@ def __init__(self, elements, labels=None, name=None, longname=None, real=False,
                 if elshape is None: elshape = el.shape
                 else: assert(elshape == el.shape), "Inconsistent element shapes!"
                 self.elements.append(el)
-            dim = int(_np.product(elshape))
+            dim = int(_np.prod(elshape))
         self.ellookup = {lbl: el for lbl, el in zip(self.labels, self.elements)}  # fast by-label element lookup
 
         if vector_elements is not None:
@@ -1459,7 +1459,7 @@ def to_elementstd_transform_matrix(self):
             number of vectors).
         """
         assert(not self.sparse), "to_elementstd_transform_matrix not implemented for sparse mode"
-        expanddim = self.elsize  # == _np.product(self.elshape)
+        expanddim = self.elsize  # == _np.prod(self.elshape)
         if self.sparse:
             toSimpleStd = _sps.lil_matrix((expanddim, self.size), dtype='complex')
         else:
@@ -1617,7 +1617,7 @@ def dim(self):
         spans.  Equivalently, the length of the `vector_elements` of the
         basis.
         """
-        dim = int(_np.product([c.dim for c in self.component_bases]))
+        dim = int(_np.prod([c.dim for c in self.component_bases]))
 
         #NOTE: this is actually to restrictive -- what we need is a test/flag for whether the elements of a
         # basis are in their "natrual" representation where it makes sense to take tensor products.  For
@@ -1635,7 +1635,7 @@ def size(self):
         """
         The number of elements (or vector-elements) in the basis.
         """
-        return int(_np.product([c.size for c in self.component_bases]))
+        return int(_np.prod([c.size for c in self.component_bases]))
 
     @property
     def elshape(self):
diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 2f4f92bf8..432b81c75 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -194,7 +194,7 @@ def _count_uptriangle_labels_for_support(cls, support, left_support, type_str, t
             right_offsets = [(i + 1 if ii < ifirst_trivial else 0) for ii, i in enumerate(left_inds)]
             if n1 == n: right_offsets[-1] += 1  # advance past diagonal element
             start_at = _np.dot(right_offsets, placevals)
-            cnt += _np.product(right_lengths) - start_at
+            cnt += _np.prod(right_lengths) - start_at
 
         return cnt
 
diff --git a/pygsti/baseobjs/opcalc/slowopcalc.py b/pygsti/baseobjs/opcalc/slowopcalc.py
index 81f55a418..db2af2142 100644
--- a/pygsti/baseobjs/opcalc/slowopcalc.py
+++ b/pygsti/baseobjs/opcalc/slowopcalc.py
@@ -259,4 +259,4 @@ def compact_deriv(vtape, ctape, wrt_params):
 
 
 def float_product(ar):
-    return _np.product(ar)
+    return _np.prod(ar)
diff --git a/pygsti/baseobjs/polynomial.py b/pygsti/baseobjs/polynomial.py
index eb77c2e16..f82c4b8b1 100644
--- a/pygsti/baseobjs/polynomial.py
+++ b/pygsti/baseobjs/polynomial.py
@@ -278,7 +278,7 @@ def evaluate(self, variable_values):
         #FUTURE: make this function smarter (Russian peasant)
         ret = 0
         for ivar, coeff in self.coeffs.items():
-            ret += coeff * _np.product([variable_values[i] for i in ivar])
+            ret += coeff * _np.prod([variable_values[i] for i in ivar])
         return ret
 
     def compact(self, complex_coeff_tape=True):
@@ -768,7 +768,7 @@ def to_rep(self):  # , max_num_vars=None not needed anymore -- given at __init__
 #         #FUTURE: make this function smarter (Russian peasant)
 #         ret = 0
 #         for ivar, coeff in self.items():
-#             ret += coeff * _np.product([variable_values[i] for i in ivar])
+#             ret += coeff * _np.prod([variable_values[i] for i in ivar])
 #         assert(_np.isclose(ret, self.fastpoly.evaluate(variable_values)))
 #         self._check_fast_polynomial()
 #         return ret
diff --git a/pygsti/baseobjs/resourceallocation.py b/pygsti/baseobjs/resourceallocation.py
index 515fb898a..db3d89489 100644
--- a/pygsti/baseobjs/resourceallocation.py
+++ b/pygsti/baseobjs/resourceallocation.py
@@ -321,7 +321,7 @@ def gather_base(self, result, local, slice_of_global, unit_ralloc=None, all_gath
                 #OLD: gathered_data = gather_comm.allgather(local)  # could change this to Allgatherv (?)
                 slices = gather_comm.allgather(slice_of_global if participating else None)
                 shapes = gather_comm.allgather(local.shape if participating else (0,))
-                sizes = [_np.product(shape) for shape in shapes]
+                sizes = [_np.prod(shape) for shape in shapes]
                 gathered_data = _np.empty(sum(sizes), dtype=local.dtype)
                 gather_comm.Allgatherv(local.flatten() if participating
                                        else _np.empty(0, dtype=local.dtype), (gathered_data, sizes))
@@ -331,7 +331,7 @@ def gather_base(self, result, local, slice_of_global, unit_ralloc=None, all_gath
                 slices = gather_comm.gather(slice_of_global if participating else None, root=0)
 
                 if gather_comm.rank == 0:
-                    sizes = [_np.product(shape) for shape in shapes]
+                    sizes = [_np.prod(shape) for shape in shapes]
                     gathered_data = _np.empty(sum(sizes), dtype=local.dtype)
                     recvbuf = (gathered_data, sizes)
                 else:
diff --git a/pygsti/baseobjs/statespace.py b/pygsti/baseobjs/statespace.py
index 91f325331..d6c3a4fe2 100644
--- a/pygsti/baseobjs/statespace.py
+++ b/pygsti/baseobjs/statespace.py
@@ -610,7 +610,7 @@ def udim(self):
         """
         Integer Hilbert (unitary operator) space dimension of this quantum state space.
         """
-        return _np.product(self.qudit_udims)
+        return _np.prod(self.qudit_udims)
 
     @property
     def dim(self):
@@ -1060,17 +1060,17 @@ def is_label(x):
         self.tpb_dims = []
         self.tpb_udims = []
         for iTPB, tpbLabels in enumerate(self.labels):
-            float_prod = _np.product(_np.array([self.label_dims[lbl] for lbl in tpbLabels], 'd'))
+            float_prod = _np.prod(_np.array([self.label_dims[lbl] for lbl in tpbLabels], 'd'))
             if float_prod >= float(_sys.maxsize):  # too many qubits to hold dimension in an integer
                 self.tpb_dims.append(_np.inf)
             else:
-                self.tpb_dims.append(int(_np.product([self.label_dims[lbl] for lbl in tpbLabels])))
+                self.tpb_dims.append(int(_np.prod([self.label_dims[lbl] for lbl in tpbLabels])))
 
-            float_prod = _np.product(_np.array([self.label_udims[lbl] for lbl in tpbLabels], 'd'))
+            float_prod = _np.prod(_np.array([self.label_udims[lbl] for lbl in tpbLabels], 'd'))
             if float_prod >= float(_sys.maxsize):  # too many qubits to hold dimension in an integer
                 self.tpb_udims.append(_np.inf)
             else:
-                self.tpb_udims.append(int(_np.product([self.label_udims[lbl] for lbl in tpbLabels])))
+                self.tpb_udims.append(int(_np.prod([self.label_udims[lbl] for lbl in tpbLabels])))
 
             self.tpb_index.update({lbl: iTPB for lbl in tpbLabels})
 
diff --git a/pygsti/evotypes/densitymx_slow/effectreps.py b/pygsti/evotypes/densitymx_slow/effectreps.py
index 7109e886d..37d0d6599 100644
--- a/pygsti/evotypes/densitymx_slow/effectreps.py
+++ b/pygsti/evotypes/densitymx_slow/effectreps.py
@@ -96,7 +96,7 @@ def __init__(self, povm_factors, effect_labels, state_space):
         self.factor_dims = factordims
         self.max_factor_dim = max_factor_dim  # Unused
         state_space = _StateSpace.cast(state_space)
-        assert(_np.product(factordims) == state_space.dim)
+        assert(_np.prod(factordims) == state_space.dim)
         super(EffectRepTensorProduct, self).__init__(state_space)
         self.factor_effects_have_changed()
 
diff --git a/pygsti/evotypes/densitymx_slow/opreps.py b/pygsti/evotypes/densitymx_slow/opreps.py
index 08d79e825..8feb14d95 100644
--- a/pygsti/evotypes/densitymx_slow/opreps.py
+++ b/pygsti/evotypes/densitymx_slow/opreps.py
@@ -329,7 +329,7 @@ def __init__(self, state_space, target_labels, embedded_rep):
         # final map just acts as identity w.r.t.
         labelIndices = [tensorProdBlkLabels.index(label) for label in target_labels]
         actionInds = _np.array(labelIndices, _np.int64)
-        assert(_np.product([numBasisEls[i] for i in actionInds]) == embedded_rep.dim), \
+        assert(_np.prod([numBasisEls[i] for i in actionInds]) == embedded_rep.dim), \
             "Embedded operation has dimension (%d) inconsistent with the given target labels (%s)" % (
                 embedded_rep.dim, str(target_labels))
 
@@ -337,7 +337,7 @@ def __init__(self, state_space, target_labels, embedded_rep):
         iActiveBlock = iTensorProdBlk
         nComponents = len(state_space.tensor_product_block_labels(iActiveBlock))
         #embeddedDim = embedded_rep.dim
-        blocksizes = _np.array([_np.product(state_space.tensor_product_block_dimensions(k))
+        blocksizes = _np.array([_np.prod(state_space.tensor_product_block_dimensions(k))
                                 for k in range(nBlocks)], _np.int64)
 
         self.embedded_rep = embedded_rep
diff --git a/pygsti/evotypes/densitymx_slow/statereps.py b/pygsti/evotypes/densitymx_slow/statereps.py
index 6bc230a2a..b5d3e4e8a 100644
--- a/pygsti/evotypes/densitymx_slow/statereps.py
+++ b/pygsti/evotypes/densitymx_slow/statereps.py
@@ -143,7 +143,7 @@ def __reduce__(self):
 class StateRepTensorProduct(StateRep):
     def __init__(self, factor_state_reps, state_space):
         self.factor_reps = factor_state_reps
-        dim = _np.product([fct.dim for fct in self.factor_reps])
+        dim = _np.prod([fct.dim for fct in self.factor_reps])
         super(StateRepTensorProduct, self).__init__(_np.zeros(dim, 'd'), state_space)
         self.reps_have_changed()
 
diff --git a/pygsti/evotypes/stabilizer_slow/statereps.py b/pygsti/evotypes/stabilizer_slow/statereps.py
index bb070fccd..b47417f08 100644
--- a/pygsti/evotypes/stabilizer_slow/statereps.py
+++ b/pygsti/evotypes/stabilizer_slow/statereps.py
@@ -104,7 +104,7 @@ class StateRepTensorProduct(StateRep):
     def __init__(self, factor_state_reps, state_space):
         self.factor_reps = factor_state_reps
         n = sum([sf.nqubits for sf in self.factor_reps])  # total number of qubits
-        np = int(_np.product([len(sf.pvectors) for sf in self.factor_reps]))
+        np = int(_np.prod([len(sf.pvectors) for sf in self.factor_reps]))
 
         super(StateRepTensorProduct, self).__init__(_np.zeros((2 * n, 2 * n), _np.int64),
                                                     _np.zeros((np, 2 * n), _np.int64),
diff --git a/pygsti/evotypes/statevec_slow/effectreps.py b/pygsti/evotypes/statevec_slow/effectreps.py
index 8c391cc7b..d1a14ebc6 100644
--- a/pygsti/evotypes/statevec_slow/effectreps.py
+++ b/pygsti/evotypes/statevec_slow/effectreps.py
@@ -98,7 +98,7 @@ def __init__(self, povm_factors, effect_labels, state_space):
         factordims = _np.ascontiguousarray(
             _np.array([fct.state_space.udim for fct in povm_factors], _np.int64))
 
-        #dim = _np.product(factordims)
+        #dim = _np.prod(factordims)
         self.povm_factors = povm_factors
         self.effect_labels = effect_labels
         self.kron_array = kron_array
@@ -106,7 +106,7 @@ def __init__(self, povm_factors, effect_labels, state_space):
         self.nfactors = len(self.povm_factors)
         self.max_factor_dim = max_factor_dim  # Unused
         state_space = _StateSpace.cast(state_space)
-        assert(_np.product(factordims) == state_space.udim)
+        assert(_np.prod(factordims) == state_space.udim)
         super(EffectRepTensorProduct, self).__init__(state_space)
         self.factor_effects_have_changed()
 
diff --git a/pygsti/evotypes/statevec_slow/opreps.py b/pygsti/evotypes/statevec_slow/opreps.py
index c66a66aa8..b60fadd51 100644
--- a/pygsti/evotypes/statevec_slow/opreps.py
+++ b/pygsti/evotypes/statevec_slow/opreps.py
@@ -203,7 +203,7 @@ def __init__(self, state_space, target_labels, embedded_rep):
         # final map just acts as identity w.r.t.
         labelIndices = [tensorProdBlkLabels.index(label) for label in target_labels]
         actionInds = _np.array(labelIndices, _np.int64)
-        assert(_np.product([numBasisEls[i] for i in actionInds]) == embedded_rep.dim), \
+        assert(_np.prod([numBasisEls[i] for i in actionInds]) == embedded_rep.dim), \
             "Embedded operation has dimension (%d) inconsistent with the given target labels (%s)" % (
                 embedded_rep.dim, str(target_labels))
 
@@ -212,7 +212,7 @@ def __init__(self, state_space, target_labels, embedded_rep):
         iActiveBlock = iTensorProdBlk
         nComponents = len(state_space.tensor_product_block_labels(iActiveBlock))
         embeddedDim = embedded_rep.dim  # a *unitary* dim - see .dim property above
-        blocksizes = _np.array([_np.product(state_space.tensor_product_block_udimensions(k))
+        blocksizes = _np.array([_np.prod(state_space.tensor_product_block_udimensions(k))
                                 for k in range(nBlocks)], _np.int64)
 
         self.target_labels = target_labels
diff --git a/pygsti/evotypes/statevec_slow/statereps.py b/pygsti/evotypes/statevec_slow/statereps.py
index b1b29bdf1..7b5682733 100644
--- a/pygsti/evotypes/statevec_slow/statereps.py
+++ b/pygsti/evotypes/statevec_slow/statereps.py
@@ -123,7 +123,7 @@ def actionable_staterep(self):
 class StateRepTensorProduct(StateRep):
     def __init__(self, factor_state_reps, state_space):
         self.factor_reps = factor_state_reps
-        dim = _np.product([fct.dim for fct in self.factor_reps])
+        dim = _np.prod([fct.dim for fct in self.factor_reps])
         # FUTURE TODO: below compute a tensorprod basis instead of punting and passing `None`
         super(StateRepTensorProduct, self).__init__(_np.zeros(dim, complex), state_space, None)
         self.reps_have_changed()
diff --git a/pygsti/extras/idletomography/pauliobjs.py b/pygsti/extras/idletomography/pauliobjs.py
index 9accc3bc0..1d6d443cc 100644
--- a/pygsti/extras/idletomography/pauliobjs.py
+++ b/pygsti/extras/idletomography/pauliobjs.py
@@ -405,7 +405,7 @@ def ri_sign(pauli1, pauli2, parity):
         sign = (-1)**((num_i + 1) / 2) * _np.prod([ri_sign(pauli1, pauli2, p)
                                                    for pauli1, pauli2, p in zip(s1, s2, parities)])
         if isinstance(other, NQPauliOp): other_sign = other.sign
-        elif isinstance(other, NQPauliState): other_sign = _np.product(other.signs)
+        elif isinstance(other, NQPauliState): other_sign = _np.prod(other.signs)
         else: raise ValueError("Can't take commutator with %s type" % str(type(other)))
 
         return NQPauliOp(op, sign * self.sign * other_sign)
diff --git a/pygsti/extras/interpygate/core.py b/pygsti/extras/interpygate/core.py
index cd0a7c232..d7b205146 100644
--- a/pygsti/extras/interpygate/core.py
+++ b/pygsti/extras/interpygate/core.py
@@ -565,7 +565,7 @@ def compute_data(self, comm=None, mpi_workers_per_process=1, verbosity=0):
 
         if rank in root_ranks:
             #Only root ranks store data (fn_to_interpolate only needs to return results on root proc)
-            flat_data = _np.empty(len(my_points) * int(_np.product(expected_fn_output_shape)), dtype='d')
+            flat_data = _np.empty(len(my_points) * int(_np.prod(expected_fn_output_shape)), dtype='d')
             data = flat_data.view(); data.shape = (len(my_points),) + expected_fn_output_shape
             if (comm is not None):
                 printer.log("Group %d processing %d points on %d processors." % (color, len(my_points),
diff --git a/pygsti/forwardsims/distforwardsim.py b/pygsti/forwardsims/distforwardsim.py
index e7c1db727..d0a32bed4 100644
--- a/pygsti/forwardsims/distforwardsim.py
+++ b/pygsti/forwardsims/distforwardsim.py
@@ -456,7 +456,7 @@ def _compute_processor_distribution(self, array_types, nprocs, num_params, num_c
             else self._pblk_sizes[0:len(param_dimensions)]  # automatically set these?
 
         if self._processor_grid is not None:
-            assert(_np.product(self._processor_grid) <= nprocs), "`processor_grid` must multiply to # of procs!"
+            assert(_np.prod(self._processor_grid) <= nprocs), "`processor_grid` must multiply to # of procs!"
             na = self._processor_grid[0]
             natoms = max(na, self._num_atoms) if (self._num_atoms is not None) else na
             npp = ()
diff --git a/pygsti/forwardsims/mapforwardsim.py b/pygsti/forwardsims/mapforwardsim.py
index 81b501f44..d86b9e9a3 100644
--- a/pygsti/forwardsims/mapforwardsim.py
+++ b/pygsti/forwardsims/mapforwardsim.py
@@ -250,9 +250,9 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
             array_types, nprocs, num_params, len(circuits), default_natoms=2 * self.model.dim)  # heuristic?
 
         printer.log("MapLayout: %d processors divided into %s (= %d) grid along circuit and parameter directions." %
-                    (nprocs, ' x '.join(map(str, (na,) + npp)), _np.product((na,) + npp)))
+                    (nprocs, ' x '.join(map(str, (na,) + npp)), _np.prod((na,) + npp)))
         printer.log("   %d atoms, parameter block size limits %s" % (natoms, str(param_blk_sizes)))
-        assert(_np.product((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
+        assert(_np.prod((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
 
         layout = _MapCOPALayout(circuits, self.model, dataset, self._max_cache_size, natoms, na, npp,
                                 param_dimensions, param_blk_sizes, resource_alloc, verbosity)
diff --git a/pygsti/forwardsims/matrixforwardsim.py b/pygsti/forwardsims/matrixforwardsim.py
index 21267a493..d82543b3e 100644
--- a/pygsti/forwardsims/matrixforwardsim.py
+++ b/pygsti/forwardsims/matrixforwardsim.py
@@ -1093,9 +1093,9 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
             na, npp = 1, (1, 1)  # save all processor division for within the (single) atom, for different timestamps
 
         printer.log("MatrixLayout: %d processors divided into %s (= %d) grid along circuit and parameter directions." %
-                    (nprocs, ' x '.join(map(str, (na,) + npp)), _np.product((na,) + npp)))
+                    (nprocs, ' x '.join(map(str, (na,) + npp)), _np.prod((na,) + npp)))
         printer.log("   %d atoms, parameter block size limits %s" % (natoms, str(param_blk_sizes)))
-        assert(_np.product((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
+        assert(_np.prod((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
 
         layout = _MatrixCOPALayout(circuits, self.model, dataset, natoms,
                                    na, npp, param_dimensions, param_blk_sizes, resource_alloc, verbosity)
diff --git a/pygsti/forwardsims/termforwardsim.py b/pygsti/forwardsims/termforwardsim.py
index ea07fe1b2..4cd23791c 100644
--- a/pygsti/forwardsims/termforwardsim.py
+++ b/pygsti/forwardsims/termforwardsim.py
@@ -326,9 +326,9 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
             array_types, nprocs, num_params, len(circuits), default_natoms=nprocs)
 
         printer.log("TermLayout: %d processors divided into %s (= %d) grid along circuit and parameter directions." %
-                    (nprocs, ' x '.join(map(str, (na,) + npp)), _np.product((na,) + npp)))
+                    (nprocs, ' x '.join(map(str, (na,) + npp)), _np.prod((na,) + npp)))
         printer.log("   %d atoms, parameter block size limits %s" % (natoms, str(param_blk_sizes)))
-        assert(_np.product((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
+        assert(_np.prod((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
 
         layout = _TermCOPALayout(circuits, self.model, dataset, natoms, na, npp, param_dimensions,
                                  param_blk_sizes, resource_alloc, printer)
@@ -950,7 +950,7 @@ def _achieved_and_max_sopm_jacobian_atom(self, layout_atom):
             Eops = [self.model.circuit_layer_operator(elbl, 'povm') for elbl in elabels]
             partial_op_maxmag_values = [op.total_term_magnitude() for op in partial_ops]
             Eop_maxmag_values = [Eop.total_term_magnitude() for Eop in Eops]
-            maxmag_partial_product = _np.product(partial_op_maxmag_values)
+            maxmag_partial_product = _np.prod(partial_op_maxmag_values)
             maxmag_products = [maxmag_partial_product * Eop_val for Eop_val in Eop_maxmag_values]
 
             deriv = _np.zeros((len(elabels), Np), 'd')
diff --git a/pygsti/forwardsims/termforwardsim_calc_generic.py b/pygsti/forwardsims/termforwardsim_calc_generic.py
index 792bf2606..c23b35d30 100644
--- a/pygsti/forwardsims/termforwardsim_calc_generic.py
+++ b/pygsti/forwardsims/termforwardsim_calc_generic.py
@@ -72,7 +72,7 @@ def prs_as_polynomials(fwdsim, rholabel, elabels, circuit, polynomial_vindices_p
     # use get_direct_order_terms(order, order_base) w/order_base=0.1(?) instead of taylor_order_terms??
     # below: replace prps with: prs = _np.zeros(len(elabels),complex)  # an array in "bulk" mode
     #  use *= or * instead of .mult( and .scale(
-    #  e.g. res = _np.product([f.coeff for f in factors])
+    #  e.g. res = _np.prod([f.coeff for f in factors])
     #       res *= (pLeft * pRight)
     # - add assert(_np.linalg.norm(_np.imag(prs)) < 1e-6) at end and return _np.real(prs)
 
@@ -227,7 +227,7 @@ def prs_as_polynomials(fwdsim, rholabel, elabels, circuit, polynomial_vindices_p
 
             # #DEBUG!!!
             # db_nfactors = [len(l) for l in factor_lists]
-            # db_totfactors = _np.product(db_nfactors)
+            # db_totfactors = _np.prod(db_nfactors)
             # db_factor_cnt += db_totfactors
             # DEBUG_FCOUNT += db_totfactors
             # db_part_cnt += 1
@@ -347,7 +347,7 @@ def circuit_achieved_and_max_sopm(fwdsim, rholabel, elabels, circuit, repcache,
 
     ops = [fwdsim.model._circuit_layer_operator(rholabel, 'prep')] + \
         [fwdsim.model._circuit_layer_operator(glbl, 'op') for glbl in circuit]
-    max_sum_of_pathmags = _np.product([op.total_term_magnitude for op in ops])
+    max_sum_of_pathmags = _np.prod([op.total_term_magnitude for op in ops])
     max_sum_of_pathmags = _np.array(
         [max_sum_of_pathmags * fwdsim.model._circuit_layer_operator(elbl, 'povm').total_term_magnitude
          for elbl in elabels], 'd')
@@ -459,7 +459,7 @@ def find_best_pathmagnitude_threshold(fwdsim, rholabel, elabels, circuit, polyno
 
     ops = [fwdsim.model._circuit_layer_operator(rholabel, 'prep')] + \
         [fwdsim.model._circuit_layer_operator(glbl, 'op') for glbl in circuit]
-    max_sum_of_pathmags = _np.product([op.total_term_magnitude for op in ops])
+    max_sum_of_pathmags = _np.prod([op.total_term_magnitude for op in ops])
     max_sum_of_pathmags = _np.array(
         [max_sum_of_pathmags * fwdsim.model._circuit_layer_operator(elbl, 'povm').total_term_magnitude
          for elbl in elabels], 'd')
@@ -838,7 +838,7 @@ def _prs_as_pruned_polys(fwdsim, rholabel, elabels, circuit, repcache, comm=None
 
     ops = [fwdsim.model._circuit_layer_operator(rholabel, 'prep')] + \
         [fwdsim.model._circuit_layer_operator(glbl, 'op') for glbl in circuit]
-    max_sum_of_pathmags = _np.product([op.total_term_magnitude for op in ops])
+    max_sum_of_pathmags = _np.prod([op.total_term_magnitude for op in ops])
     max_sum_of_pathmags = _np.array(
         [max_sum_of_pathmags * fwdsim.model._circuit_layer_operator(elbl, 'povm').total_term_magnitude
          for elbl in elabels], 'd')
diff --git a/pygsti/modelmembers/operations/composedop.py b/pygsti/modelmembers/operations/composedop.py
index 4e2574d74..c8e941cab 100644
--- a/pygsti/modelmembers/operations/composedop.py
+++ b/pygsti/modelmembers/operations/composedop.py
@@ -543,7 +543,7 @@ def taylor_order_terms_above_mag(self, order, max_polynomial_vars, min_term_mag)
             #                 for i, pi in enumerate(p)]
             factor_lists = [factor_lists_cache[i][pi] for i, pi in enumerate(p)]
             for factors in _itertools.product(*factor_lists):
-                mag = _np.product([factor.magnitude for factor in factors])
+                mag = _np.prod([factor.magnitude for factor in factors])
                 if mag >= min_term_mag:
                     terms.append(_term.compose_terms_with_mag(factors, mag))
         return terms
@@ -585,7 +585,7 @@ def total_term_magnitude(self):
         #  of an errorgen or operator.
         # In this case, since the taylor expansions are composed (~multiplied),
         # the total term magnitude is just the product of those of the components.
-        return _np.product([f.total_term_magnitude for f in self.factorops])
+        return _np.prod([f.total_term_magnitude for f in self.factorops])
 
     @property
     def total_term_magnitude_deriv(self):
@@ -601,7 +601,7 @@ def total_term_magnitude_deriv(self):
             An array of length self.num_params
         """
         opmags = [f.total_term_magnitude for f in self.factorops]
-        product = _np.product(opmags)
+        product = _np.prod(opmags)
         ret = _np.zeros(self.num_params, 'd')
         for opmag, f, f_local_inds in zip(opmags, self.factorops, self._submember_rpindices):
             #f_local_inds = _modelmember._decompose_gpindices(
diff --git a/pygsti/modelmembers/operations/embeddedop.py b/pygsti/modelmembers/operations/embeddedop.py
index 3e430a94d..c21ce8397 100644
--- a/pygsti/modelmembers/operations/embeddedop.py
+++ b/pygsti/modelmembers/operations/embeddedop.py
@@ -162,9 +162,9 @@ def _iter_matrix_elements_precalc(self, on_space):
 
         # number of basis elements preceding our block's elements
         if on_space == "Hilbert":
-            blockDims = [_np.product(tpb_dims) for tpb_dims in self.state_space.tensor_product_blocks_udimensions]
+            blockDims = [_np.prod(tpb_dims) for tpb_dims in self.state_space.tensor_product_blocks_udimensions]
         else:
-            blockDims = [_np.product(tpb_dims) for tpb_dims in self.state_space.tensor_product_blocks_dimensions]
+            blockDims = [_np.prod(tpb_dims) for tpb_dims in self.state_space.tensor_product_blocks_dimensions]
         offset = sum(blockDims[0:iTensorProdBlk])
 
         return divisors, multipliers, sorted_bili, basisInds_noop, offset
diff --git a/pygsti/modelmembers/povms/tensorprodeffect.py b/pygsti/modelmembers/povms/tensorprodeffect.py
index afa765753..57c59da65 100644
--- a/pygsti/modelmembers/povms/tensorprodeffect.py
+++ b/pygsti/modelmembers/povms/tensorprodeffect.py
@@ -337,7 +337,7 @@ def deriv_wrt_params(self, wrt_filter=None):
         #HACK to deal with fact that output of to_dense is really what is differentiated
         # but this may not match self.dim == self.state_space.dim, e.g. for pure state vecs.
         dims = [len(fct.to_dense(on_space='minimal')) for fct in self.factors]
-        dim = int(_np.product(dims))
+        dim = int(_np.prod(dims))
 
         derivMx = _np.zeros((dim, self.num_params), typ)
 
diff --git a/pygsti/modelmembers/povms/tensorprodpovm.py b/pygsti/modelmembers/povms/tensorprodpovm.py
index 7750c5874..20c008e52 100644
--- a/pygsti/modelmembers/povms/tensorprodpovm.py
+++ b/pygsti/modelmembers/povms/tensorprodpovm.py
@@ -43,7 +43,7 @@ class TensorProductPOVM(_POVM):
     """
 
     def __init__(self, factor_povms, evotype="auto", state_space=None):
-        dim = _np.product([povm.state_space.dim for povm in factor_povms])
+        dim = _np.prod([povm.state_space.dim for povm in factor_povms])
         if state_space is None:
             state_space = _statespace.default_space_for_dim(dim)
         else:
@@ -102,7 +102,7 @@ def __iter__(self):
         return self.keys()
 
     def __len__(self):
-        return _np.product([len(fk) for fk in self._factor_keys])
+        return _np.prod([len(fk) for fk in self._factor_keys])
 
     def keys(self):
         """
diff --git a/pygsti/modelmembers/states/tensorprodstate.py b/pygsti/modelmembers/states/tensorprodstate.py
index 8583e4f15..c71e59630 100644
--- a/pygsti/modelmembers/states/tensorprodstate.py
+++ b/pygsti/modelmembers/states/tensorprodstate.py
@@ -285,7 +285,7 @@ def deriv_wrt_params(self, wrt_filter=None):
         #HACK to deal with fact that output of to_dense is really what is differentiated
         # but this may not match self.dim == self.state_space.dim, e.g. for pure state vecs.
         dims = [len(fct.to_dense(on_space='minimal')) for fct in self.factors]
-        dim = int(_np.product(dims))
+        dim = int(_np.prod(dims))
 
         derivMx = _np.zeros((dim, self.num_params), typ)
 
diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py
index cd620e9e8..e2a311a29 100644
--- a/pygsti/models/explicitmodel.py
+++ b/pygsti/models/explicitmodel.py
@@ -1504,7 +1504,7 @@ def extract_unitary(Umx, U_sslbls, extracted_sslbls):
             # assume this is a kronecker product (check this in FUTURE?), so just fill extracted
             # unitary by fixing all non-extracted qudits (assumed identity-action on these) to 0
             # and looping over extracted ones:
-            U_extracted = _np.zeros((_np.product(extracted_udims), _np.product(extracted_udims)), complex)
+            U_extracted = _np.zeros((_np.prod(extracted_udims), _np.prod(extracted_udims)), complex)
             for ii, itup in enumerate(_itertools.product(*[range(ud) for ud in extracted_udims])):
                 i = _np.dot(extracted_inc, itup)
                 for jj, jtup in enumerate(_itertools.product(*[range(ud) for ud in extracted_udims])):
diff --git a/pygsti/models/modelconstruction.py b/pygsti/models/modelconstruction.py
index dc53a906c..1ef1ba189 100644
--- a/pygsti/models/modelconstruction.py
+++ b/pygsti/models/modelconstruction.py
@@ -240,7 +240,7 @@ def to_labels(lbls):
         if opName == "I":
             # qubit labels (TODO: what about 'L' labels? -- not sure if they work with this...)
             labels = to_labels(args)
-            stateSpaceUDim = int(_np.product([state_space.label_udimension(l) for l in labels]))
+            stateSpaceUDim = int(_np.prod([state_space.label_udimension(l) for l in labels]))
             # a complex 2x2 mx unitary for the identity in Pauli-product basis
             Uop = _op.StaticUnitaryOp(_np.identity(stateSpaceUDim, 'complex'), 'pp', build_evotype)
 
diff --git a/pygsti/serialization/jsoncodec.py b/pygsti/serialization/jsoncodec.py
index aa7594c22..1f9cd18a2 100644
--- a/pygsti/serialization/jsoncodec.py
+++ b/pygsti/serialization/jsoncodec.py
@@ -231,7 +231,7 @@ def _encode_std_obj(py_obj, binary):
         if kind == 'O':
             #Special case of object arrays:  store flattened array data
             data = [encode_obj(el, binary) for el in py_obj.flat]
-            assert(len(data) == _np.product(py_obj.shape))
+            assert(len(data) == _np.prod(py_obj.shape))
         else:
             data = py_obj.tobytes() if binary else _tostr(_base64.b64encode(py_obj.tobytes()))
 
diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index 6894adab6..bbc8bed1c 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -103,8 +103,8 @@ def _embed(mx, target_labels, state_space):  # BOTTLENECK
             gen2_expanded = _embed(gen2, gen2_sslbls, action_space)  # embed gen2 into action_space
             gen2_expanded *= scale  # so that gen2_expanded is an embedded *dual* generator
             if _sps.issparse(gen2_expanded):
-                flat_gen2_expanded = gen2_expanded.reshape((1, _np.product(gen2_expanded.shape)))
-                flat_gauge_action_deriv = gauge_action_deriv.reshape((_np.product(gauge_action_deriv.shape), 1))
+                flat_gen2_expanded = gen2_expanded.reshape((1, _np.prod(gen2_expanded.shape)))
+                flat_gauge_action_deriv = gauge_action_deriv.reshape((_np.prod(gauge_action_deriv.shape), 1))
                 val = flat_gen2_expanded.dot(flat_gauge_action_deriv)[0, 0]  # Note: gen2 is a *dual* generator
             else:
                 val = _np.vdot(gen2_expanded.flat, gauge_action_deriv.flat)
diff --git a/pygsti/tools/mpitools.py b/pygsti/tools/mpitools.py
index d47e119fc..36bcac634 100644
--- a/pygsti/tools/mpitools.py
+++ b/pygsti/tools/mpitools.py
@@ -1034,8 +1034,8 @@ def processor_group_size(nprocs, number_of_tasks):
         fctrs = sorted(_prime_factors(nprocs)); i = 1
         if int(_np.ceil(number_of_tasks)) in fctrs:
             return int(_np.ceil(number_of_tasks))  # we got lucky
-        while _np.product(fctrs[0:i]) < number_of_tasks: i += 1
-        return _np.product(fctrs[0:i])
+        while _np.prod(fctrs[0:i]) < number_of_tasks: i += 1
+        return _np.prod(fctrs[0:i])
 
 
 def sum_arrays(local_array, owners, comm):
diff --git a/pygsti/tools/sharedmemtools.py b/pygsti/tools/sharedmemtools.py
index f5f956c03..8c81b10ec 100644
--- a/pygsti/tools/sharedmemtools.py
+++ b/pygsti/tools/sharedmemtools.py
@@ -129,7 +129,7 @@ def create_shared_ndarray(resource_alloc, shape, dtype, zero_out=False, memory_t
         :function:`cleanup_shared_ndarray` to ensure `ar` is deallocated properly.
     """
     hostcomm = resource_alloc.host_comm if shared_mem_is_enabled() else None
-    nelements = _np.product(shape)
+    nelements = _np.prod(shape)
     if hostcomm is None or nelements == 0:  # Note: shared memory must be for size > 0
         # every processor allocates its own memory
         if memory_tracker is not None: memory_tracker.add_tracked_memory(nelements)