Merge branch 'master' of https://github.com/quantumlib/cirq into v0.1…

…0.0-dev

.github/CODEOWNERS

@@ -8,4 +8,5 @@
 
 **/pasqal/**/*.* @HGSilveri @quantumlib/cirq-maintainers @vtomole @cduck
 
-cirq/experiments/**/*.* @mrwojtek @quantumlib/cirq-maintainers @vtomole @cduck
+cirq/experiments/**/*.* @mrwojtek @quantumlib/cirq-maintainers @vtomole @cduck
+cirq/docs/qcvv/**/*.* @mrwojtek @quantumlib/cirq-maintainers @vtomole @cduck

cirq/experiments/xeb_fitting.py

@@ -206,6 +206,14 @@ def get_initial_simplex_and_names(
 
         return initial_simplex, names
 
+    def get_parameterized_gate(self):
+        theta = THETA_SYMBOL if self.characterize_theta else self.theta_default
+        zeta = ZETA_SYMBOL if self.characterize_zeta else self.zeta_default
+        chi = CHI_SYMBOL if self.characterize_chi else self.chi_default
+        gamma = GAMMA_SYMBOL if self.characterize_gamma else self.gamma_default
+        phi = PHI_SYMBOL if self.characterize_phi else self.phi_default
+        return ops.PhasedFSimGate(theta=theta, zeta=zeta, chi=chi, gamma=gamma, phi=phi)
+
 
 @dataclass(frozen=True)
 class SqrtISwapXEBOptions(XEBPhasedFSimCharacterizationOptions):
@@ -221,14 +229,6 @@ class SqrtISwapXEBOptions(XEBPhasedFSimCharacterizationOptions):
     def should_parameterize(op: 'cirq.Operation') -> bool:
         return op.gate == SQRT_ISWAP
 
-    def get_parameterized_gate(self):
-        theta = THETA_SYMBOL if self.characterize_theta else self.theta_default
-        zeta = ZETA_SYMBOL if self.characterize_zeta else self.zeta_default
-        chi = CHI_SYMBOL if self.characterize_chi else self.chi_default
-        gamma = GAMMA_SYMBOL if self.characterize_gamma else self.gamma_default
-        phi = PHI_SYMBOL if self.characterize_phi else self.phi_default
-        return ops.PhasedFSimGate(theta=theta, zeta=zeta, chi=chi, gamma=gamma, phi=phi)
-
 
 def parameterize_circuit(
     circuit: 'cirq.Circuit',
@@ -378,6 +378,10 @@ def characterize_phased_fsim_parameters_with_xeb_by_pair(
 
     This is appropriate if you have run parallel XEB on multiple pairs of qubits.
 
+    The optimization is done per-pair. If you have the same pair in e.g. two different
+    layers the characterization optimization will lump the data together. This is in contrast
+    with the benchmarking functionality, which will always index on `(layer_i, pair_i, pair)`.
+
     Args:
         sampled_df: The DataFrame of sampled two-qubit probability distributions returned
             from `sample_2q_xeb_circuits`.
@@ -438,7 +442,9 @@ def exponential_decay(cycle_depths: np.ndarray, a: float, layer_fid: float) -> n
     return a * layer_fid ** cycle_depths
 
 
-def _fit_exponential_decay(cycle_depths: np.ndarray, fidelities: np.ndarray) -> Tuple[float, float]:
+def _fit_exponential_decay(
+    cycle_depths: np.ndarray, fidelities: np.ndarray
+) -> Tuple[float, float, float, float]:
     """Fit an exponential model fidelity = a * layer_fid**x using nonlinear least squares.
 
     This uses `exponential_decay` as the function to fit with parameters `a` and `layer_fid`.
@@ -453,22 +459,40 @@ def _fit_exponential_decay(cycle_depths: np.ndarray, fidelities: np.ndarray) ->
         a: The first fit parameter that scales the exponential function, perhaps accounting for
             state prep and measurement (SPAM) error.
         layer_fid: The second fit parameters which serves as the base of the exponential.
+        a_std: The standard deviation of the `a` parameter estimate.
+        layer_fid_std: The standard deviation of the `layer_fid` parameter estimate.
     """
     cycle_depths = np.asarray(cycle_depths)
     fidelities = np.asarray(fidelities)
 
-    # Get initial guess by linear least squares with logarithm of model
+    # Get initial guess by linear least squares with logarithm of model.
+    # This only works for positive fidelities. We use numpy fancy indexing
+    # with `positives` (an ndarray of bools).
     positives = fidelities > 0
+    if np.sum(positives) <= 1:
+        # The sum of the boolean array is the number of `True` entries.
+        # For one or fewer positive values, we cannot perform the linear fit.
+        return 0, 0, np.inf, np.inf
     cycle_depths_pos = cycle_depths[positives]
     log_fidelities = np.log(fidelities[positives])
     slope, intercept, _, _, _ = scipy.stats.linregress(cycle_depths_pos, log_fidelities)
     layer_fid_0 = np.clip(np.exp(slope), 0, 1)
     a_0 = np.clip(np.exp(intercept), 0, 1)
 
-    (a, layer_fid), _ = scipy.optimize.curve_fit(
-        exponential_decay, cycle_depths, fidelities, p0=(a_0, layer_fid_0), bounds=((0, 0), (1, 1))
-    )
-    return a, layer_fid
+    try:
+        (a, layer_fid), pcov = scipy.optimize.curve_fit(
+            exponential_decay,
+            cycle_depths,
+            fidelities,
+            p0=(a_0, layer_fid_0),
+            bounds=((0, 0), (1, 1)),
+        )
+    except ValueError:  # coverage: ignore
+        # coverage: ignore
+        return 0, 0, np.inf, np.inf
+
+    a_std, layer_fid_std = np.sqrt(np.diag(pcov))
+    return a, layer_fid, a_std, layer_fid_std
 
 
 def _one_unique(df, name, default):
@@ -494,21 +518,26 @@ def fit_exponential_decays(fidelities_df: pd.DataFrame) -> pd.DataFrame:
         for the fit parameters "a" and "layer_fid"; and nested "cycles_depths" and "fidelities"
         lists (now grouped by pair).
     """
-    records = []
-    for pair in fidelities_df['pair'].unique():
-        f1 = fidelities_df[fidelities_df['pair'] == pair]
-        a, layer_fid = _fit_exponential_decay(f1['cycle_depth'], f1['fidelity'])
+
+    def _per_pair(f1):
+        a, layer_fid, a_std, layer_fid_std = _fit_exponential_decay(
+            f1['cycle_depth'], f1['fidelity']
+        )
         record = {
-            'pair': pair,
             'a': a,
             'layer_fid': layer_fid,
             'cycle_depths': f1['cycle_depth'].values,
             'fidelities': f1['fidelity'].values,
-            'layer_i': _one_unique(f1, 'layer_i', default=0),
-            'pair_i': _one_unique(f1, 'pair_i', default=0),
+            'a_std': a_std,
+            'layer_fid_std': layer_fid_std,
         }
-        records.append(record)
-    return pd.DataFrame(records).set_index(['pair', 'layer_i', 'pair_i'])
+        return pd.Series(record)
+
+    if 'layer_i' in fidelities_df.columns:
+        groupby = ['layer_i', 'pair_i', 'pair']
+    else:
+        groupby = ['pair']
+    return fidelities_df.groupby(groupby).apply(_per_pair)
 
 
 def before_and_after_characterization(
@@ -531,13 +560,19 @@ def before_and_after_characterization(
     fit_decay_df_c = fit_exponential_decays(characterization_result.fidelities_df)
 
     joined_df = fit_decay_df_0.join(fit_decay_df_c, how='outer', lsuffix='_0', rsuffix='_c')
+    # Remove (layer_i, pair_i) from the index. While we keep this for `fit_exponential_decays`
+    # so the same pair can be benchmarked in different contexts, the multi-pair characterization
+    # function only keys on the pair identity. This can be seen acutely by the
+    # `characterization_result.final_params` dictionary being keyed only by the pair.
+    joined_df = joined_df.reset_index().set_index('pair')
+
     joined_df['characterized_angles'] = [
-        characterization_result.final_params[pair] for pair, _, _ in joined_df.index
+        characterization_result.final_params[pair] for pair in joined_df.index
     ]
     # Take any `final_params` (for any pair). We just need the angle names.
     fp, *_ = characterization_result.final_params.values()
     for angle_name in fp.keys():
         joined_df[angle_name] = [
-            characterization_result.final_params[pair][angle_name] for pair, _, _ in joined_df.index
+            characterization_result.final_params[pair][angle_name] for pair in joined_df.index
         ]
     return joined_df

cirq/experiments/xeb_fitting_test.py

@@ -300,5 +300,19 @@ def test_fit_exponential_decays():
     rs = np.random.RandomState(999)
     cycle_depths = np.arange(3, 100, 11)
     fidelities = 0.95 * 0.98 ** cycle_depths + rs.normal(0, 0.2)
-    a, layer_fid = _fit_exponential_decay(cycle_depths, fidelities)
+    a, layer_fid, a_std, layer_fid_std = _fit_exponential_decay(cycle_depths, fidelities)
     np.testing.assert_allclose([a, layer_fid], [0.95, 0.98], atol=0.02)
+    assert 0 < a_std < 0.2 / len(cycle_depths)
+    assert 0 < layer_fid_std < 1e-3
+
+
+def test_fit_exponential_decays_negative_fids():
+    rs = np.random.RandomState(999)
+    cycle_depths = np.arange(3, 100, 11)
+    fidelities = 0.5 * 0.5 ** cycle_depths + rs.normal(0, 0.2) - 0.5
+    assert np.sum(fidelities > 0) <= 1, 'they go negative'
+    a, layer_fid, a_std, layer_fid_std = _fit_exponential_decay(cycle_depths, fidelities)
+    assert a == 0
+    assert layer_fid == 0
+    assert a_std == np.inf
+    assert layer_fid_std == np.inf

cirq/experiments/xeb_sampling.py

@@ -13,6 +13,9 @@
 # limitations under the License.
 """Estimation of fidelity associated with experimental circuit executions."""
 import concurrent
+import os
+import time
+import uuid
 from concurrent.futures.thread import ThreadPoolExecutor
 from dataclasses import dataclass
 from typing import (
@@ -32,7 +35,7 @@
 import pandas as pd
 import tqdm
 
-from cirq import ops, devices
+from cirq import ops, devices, value, protocols
 from cirq.circuits import Circuit
 from cirq.experiments.random_quantum_circuit_generation import CircuitLibraryCombination
 
@@ -78,6 +81,7 @@ def __init__(
     def __call__(self, tasks: List[_Sample2qXEBTask]) -> List[Dict[str, Any]]:
         prepared_circuits = [task.prepared_circuit for task in tasks]
         results = self.sampler.run_batch(prepared_circuits, repetitions=self.repetitions)
+        timestamp = time.time()
         assert len(results) == len(tasks)
         records = []
         for task, nested_result in zip(tasks, results):
@@ -93,6 +97,7 @@ def __call__(self, tasks: List[_Sample2qXEBTask]) -> List[Dict[str, Any]]:
                         'circuit_i': circuit_i,
                         'cycle_depth': task.cycle_depth,
                         'sampled_probs': sampled_probs,
+                        'timestamp': timestamp,
                         # Additional metadata to track *how* this circuit
                         # was zipped and executed.
                         'layer_i': task.layer_i,
@@ -266,6 +271,7 @@ def _execute_sample_2q_xeb_tasks_in_batches(
     repetitions: int,
     batch_size: int,
     progress_bar: Callable[..., ContextManager],
+    dataset_directory: Optional[str] = None,
 ) -> List[Dict[str, Any]]:
     """Helper function used in `sample_2q_xeb_circuits` to batch and execute sampling tasks."""
     n_tasks = len(tasks)
@@ -278,9 +284,13 @@ def _execute_sample_2q_xeb_tasks_in_batches(
         futures = [pool.submit(run_batch, task_batch) for task_batch in batched_tasks]
 
         records = []
-        with progress_bar(total=n_tasks) as progress:
+        with progress_bar(total=len(batched_tasks) * batch_size) as progress:
             for future in concurrent.futures.as_completed(futures):
-                records += future.result()
+                new_records = future.result()
+                if dataset_directory is not None:
+                    os.makedirs(f'{dataset_directory}', exist_ok=True)
+                    protocols.to_json(new_records, f'{dataset_directory}/xeb.{uuid.uuid4()}.json')
+                records.extend(new_records)
                 progress.update(batch_size)
     return records
 
@@ -294,6 +304,8 @@ def sample_2q_xeb_circuits(
     batch_size: int = 9,
     progress_bar: Optional[Callable[..., ContextManager]] = tqdm.tqdm,
     combinations_by_layer: Optional[List[CircuitLibraryCombination]] = None,
+    shuffle: Optional['cirq.RANDOM_STATE_OR_SEED_LIKE'] = None,
+    dataset_directory: Optional[str] = None,
 ):
     """Sample two-qubit XEB circuits given a sampler.
 
@@ -314,6 +326,11 @@ def sample_2q_xeb_circuits(
             by `circuits` will be sampled verbatim, resulting in isolated XEB characterization.
             Otherwise, this contains all the random combinations and metadata required to combine
             the circuits in `circuits` into wide, parallel-XEB-style circuits for execution.
+        shuffle: If provided, use this random state or seed to shuffle the order in which tasks
+            are executed.
+        dataset_directory: If provided, save each batch of sampled results to a file
+            `{dataset_directory}/xeb.{uuid4()}.json` where uuid4() is a random string. This can be
+            used to incrementally save results to be analyzed later.
 
     Returns:
         A pandas dataframe with index given by ['circuit_i', 'cycle_depth'].
@@ -338,6 +355,9 @@ def sample_2q_xeb_circuits(
 
     # Construct truncated-with-measurement circuits to run.
     tasks = _generate_sample_2q_xeb_tasks(zipped_circuits, cycle_depths)
+    if shuffle is not None:
+        shuffle = value.parse_random_state(shuffle)
+        shuffle.shuffle(tasks)
 
     # Batch and run tasks.
     records = _execute_sample_2q_xeb_tasks_in_batches(
@@ -347,6 +367,7 @@ def sample_2q_xeb_circuits(
         repetitions=repetitions,
         batch_size=batch_size,
         progress_bar=progress_bar,
+        dataset_directory=dataset_directory,
     )
 
     # Set up the dataframe.

cirq/experiments/xeb_sampling_test.py

@@ -11,11 +11,13 @@
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
+import glob
 import itertools
 from typing import Iterable
 
 import networkx as nx
 import numpy as np
+import pandas as pd
 import pytest
 
 import cirq
@@ -43,6 +45,7 @@ def test_sample_2q_xeb_circuits():
         sampler=cirq.Simulator(),
         circuits=circuits,
         cycle_depths=cycle_depths,
+        shuffle=np.random.RandomState(10),
     )
     assert len(df) == len(cycle_depths) * len(circuits)
     for (circuit_i, cycle_depth), row in df.iterrows():
@@ -89,11 +92,24 @@ def _manhattan_distance(qubit1: cirq.GridQubit, qubit2: cirq.GridQubit) -> int:
     )
 
 
-def test_sample_2q_parallel_xeb_circuits():
+def _assert_frame_approx_equal(df, df2, *, atol):
+    assert len(df) == len(df2)
+    for (i1, row1), (i2, row2) in zip(df.sort_index().iterrows(), df2.sort_index().iterrows()):
+        assert i1 == i2
+        for k in set(row1.keys()) | set(row2.keys()):
+            v1 = row1[k]
+            v2 = row2[k]
+            if isinstance(v1, (float, np.ndarray)):
+                np.testing.assert_allclose(v1, v2, atol=atol)
+            else:
+                assert v1 == v2, k
+
+
+def test_sample_2q_parallel_xeb_circuits(tmpdir):
     circuits = rqcg.generate_library_of_2q_circuits(
         n_library_circuits=5, two_qubit_gate=cirq.ISWAP ** 0.5, max_cycle_depth=10
     )
-    cycle_depths = [10]
+    cycle_depths = [5, 10]
     graph = _gridqubits_to_graph_device(cirq.GridQubit.rect(3, 2))
     combs = rqcg.get_random_combinations_for_device(
         n_library_circuits=len(circuits),
@@ -107,7 +123,9 @@ def test_sample_2q_parallel_xeb_circuits():
         circuits=circuits,
         cycle_depths=cycle_depths,
         combinations_by_layer=combs,
+        dataset_directory=f'{tmpdir}/my_dataset',
     )
+
     n_pairs = sum(len(c.pairs) for c in combs)
     assert len(df) == len(cycle_depths) * len(circuits) * n_pairs
     for (circuit_i, cycle_depth), row in df.iterrows():
@@ -119,6 +137,17 @@ def test_sample_2q_parallel_xeb_circuits():
         assert 0 <= row['pair_i'] < 2  # in 3x2 graph, there's a max of 2 pairs per layer
     assert len(df['pair'].unique()) == 7  # seven pairs in 3x2 graph
 
+    # Test loading from dataset
+    chunks = [record for fn in glob.glob(f'{tmpdir}/my_dataset/*') for record in cirq.read_json(fn)]
+    df2 = pd.DataFrame(chunks).set_index(['circuit_i', 'cycle_depth'])
+    df2['pair'] = [tuple(row['pair']) for _, row in df2.iterrows()]
+    actual_index_names = ['layer_i', 'pair_i', 'combination_i', 'cycle_depth']
+    _assert_frame_approx_equal(
+        df.reset_index().set_index(actual_index_names),
+        df2.reset_index().set_index(actual_index_names),
+        atol=1e-5,
+    )
+
 
 def test_sample_2q_parallel_xeb_circuits_bad_circuit_library():
     circuits = rqcg.generate_library_of_2q_circuits(

cirq/google/engine/calibration.py

@@ -245,9 +245,9 @@ def heatmap(self, key: str) -> vis.Heatmap:
             )
         value_map = {self.key_to_qubits(k): self.value_to_float(v) for k, v in metrics.items()}
         if all(len(k) == 1 for k in value_map.keys()):
-            return vis.Heatmap(value_map)
+            return vis.Heatmap(value_map, title=key.replace('_', ' ').title())
         elif all(len(k) == 2 for k in value_map.keys()):
-            return vis.TwoQubitInteractionHeatmap(value_map)
+            return vis.TwoQubitInteractionHeatmap(value_map, title=key.replace('_', ' ').title())
         raise ValueError(
             'Heatmaps are only supported if all the targets in a metric are one or two qubits.'
             + f'{key} has target qubits {value_map.keys()}'