Add a product method to MapForwardSimulator

For the generation of process matrices of circuits which is useful for germ selection and some of the reporting code (and probably other things). Also gives a warning though to discourage its use for nefarious (read: inefficient) purposes of general forward simulation.

pygsti/forwardsims/mapforwardsim.py

@@ -11,6 +11,7 @@
 #***************************************************************************************************
 
 import importlib as _importlib
+import warnings as _warnings
 
 import numpy as _np
 
@@ -639,3 +640,68 @@ def loglpp(array_to_fill, dest_indices, num_tot_outcomes, layout_atom, dataset_r
 
         return self._bulk_fill_timedep_deriv(layout, dataset, ds_circuits, num_total_outcomes,
                                              array_to_fill, dloglpp, logl_array_to_fill, loglpp)
+
+
+    #Utility method for generating process matrices for circuits. Should not be used for forward
+    #simulation when using the MapForwardSimulator.
+    def product(self, circuit, scale=False):
+        """
+        Compute the product of a specified sequence of operation labels.
+
+        Note: LinearOperator matrices are multiplied in the reversed order of the tuple. That is,
+        the first element of circuit can be thought of as the first gate operation
+        performed, which is on the far right of the product of matrices.
+
+        Parameters
+        ----------
+        circuit : Circuit or tuple of operation labels
+            The sequence of operation labels.
+
+        scale : bool, optional
+            When True, return a scaling factor (see below).
+
+        Returns
+        -------
+        product : numpy array
+            The product or scaled product of the operation matrices.
+        scale : float
+            Only returned when scale == True, in which case the
+            actual product == product * scale.  The purpose of this
+            is to allow a trace or other linear operation to be done
+            prior to the scaling.
+        """
+        _warnings.warn('Generating dense process matrix representations of circuits or gates \n'
+                       'can be inefficient and should be avoided for the purposes of forward \n'
+                       'simulation/calculation of circuit outcome probability distributions \n' 
+                       'when using the MapForwardSimulator.')
+
+        if scale:
+            scaledGatesAndExps = {}
+            scale_exp = 0
+            G = _np.identity(self.model.evotype.minimal_dim(self.model.state_space))
+            for lOp in circuit:
+                if lOp not in scaledGatesAndExps:
+                    opmx = self.model.circuit_layer_operator(lOp, 'op').to_dense(on_space='minimal')
+                    ng = max(_nla.norm(opmx), 1.0)
+                    scaledGatesAndExps[lOp] = (opmx / ng, _np.log(ng))
+
+                gate, ex = scaledGatesAndExps[lOp]
+                H = _np.dot(gate, G)   # product of gates, starting with identity
+                scale_exp += ex   # scale and keep track of exponent
+                if H.max() < _PSMALL and H.min() > -_PSMALL:
+                    nG = max(_nla.norm(G), _np.exp(-scale_exp))
+                    G = _np.dot(gate, G / nG); scale_exp += _np.log(nG)  # LEXICOGRAPHICAL VS MATRIX ORDER
+                else: G = H
+
+            old_err = _np.seterr(over='ignore')
+            scale = _np.exp(scale_exp)
+            _np.seterr(**old_err)
+
+            return G, scale
+
+        else:
+            G = _np.identity(self.model.evotype.minimal_dim(self.model.state_space))
+            for lOp in circuit:
+                G = _np.dot(self.model.circuit_layer_operator(lOp, 'op').to_dense(on_space='minimal'), G)
+                # above line: LEXICOGRAPHICAL VS MATRIX ORDER
+            return G