[WIP] Migration guide for providers to implement V2 primitives

docs/migration-guides/external-providers-primitives-v2.mdx

@@ -0,0 +1,233 @@
+---
+title: TODO
+description: Migrate to the primitives interface in external providers
+
+---
+
+# Motivation
+
+As quantum processing units (QPUs) scale and quantum algorithms evolve, 
+the Qiskit Primitives interface was introduced to provide a higher-level abstraction for 
+accessing QPUs. This interface is optimized for two core tasks in quantum algorithm develpment: 
+circuit sampling (Sampler) and expectation value estimation (Estimator).
+
+The Primitives interface has replaced the legacy ``backend.run`` model for user access of IBM Quantum Hardware,
+resulting in its deprecation and removal from the Qiskit Runtime service. A series of migration guides have
+been developed for users to upgrade their workflows to use the different primitive implementations available 
+in the qiskit ecosystem, including the runtime primitives for real hardware access, the statevector primitives
+for local statevector simulation, or the aer primitives for access to local noisy simulation.
+
+This migration guide, however, has been designed to help external providers implement their own flavor of the 
+primitives V2 interface (Sampler or Estimator) for easier pluggability into the standard qiskit workflow. These 
+implementations can be used to wrap any custom provider hardware access function 
+(for example: ``execute_workload(QPU)`` or ``resource.access()``)
+as long as the final inputs and outputs conform to the established contract from the interface.
+
+## If the external provider already implemented ``backend.run``
+
+The qiskit SDK offers ready-to-go wrappers for `backend.run` that can be easily adapted to a custom primitives workflow
+through subclassing, these are the `BackendSamplerV2` (API ref link) and `BackendEstimatorV2` (API ref link). The inputs 
+to the primitives should follow the pub syntax specified in the V2 primitives interface (link to pub docs). 
+
+An advantage of this strategy is that the wrapper can handle the input/output manipulation, so knowledge of the PUB data
+model is not required. However, this might result in a sub-optimal runtime, which could be refined through a fully
+custom primitives implementation.
+
+The snippets below show how to create a custom Estimator instance following the strategy described above. 
+The process would be analogous for a custom Sampler, modifying the base class to BackendSamplerV2.
+
+``` python
+class CustomEstimator(BackendEstimatorV2):
+    """:class:`BaseEstimatorV2 <qiskit.primitives.BaseEstimatorV2>` primitive for custom provider."""
+
+    _backend: CustomProviderResource
+
+    def __init__(
+        self,
+        backend: CustomProviderResource,
+        options: dict | None = None,
+        extra_flag_used_in_estimation: bool = True,
+        another_extra_flag: bool = False,
+    ) -> None:
+        """
+        Args:
+            backend: Custom provider resource to evaluate circuits on.
+            options: options passed to through to the underlying BackendEstimatorV2.
+            extra_flag_used_in_estimation:  if `False`, do this.
+            another_extra_flag: if `True`, do that,
+        """
+
+        # preprocess arguments if necessary according to custom flags
+        processed_backend = ...
+        processed_options = ...
+
+        super().__init__(
+            processed_backend,
+            options=processed_options,
+        )
+
+    @property
+    def backend(self) -> CustomProviderResource:
+        """Computing resource used for circuit evaluation."""
+        return self._backend
+```
+
+## If the external provider didn't implement ``backend.run`` or you look for a fully custom implementation
+
+If a new provider is developed that doesn't conform to the legacy `backend.run` interface, the pre-packaged
+wrapper may not be the most optimal route for implementing the primtives. Instead, you should implement a particular
+instance of the abstract base primitive interfaces (link to BaseSamplerV2 and BaseEstimatorV2). This route requires an
+understanding of the PUB data model for input/output handling.
+
+The following snippet shows a minimal example of an implementation of a custom Sampler primitive following this strategy.
+This example has been extracted and generalized from the StatevectorSampler implementation. The full original 
+implementation can be found in the StatevectorSampler source code.
+
+``` python
+from qiskit.primitives.base import BaseSamplerV2
+from qiskit.primitives.containers import (
+    BitArray,
+    DataBin,
+    PrimitiveResult,
+    SamplerPubResult,
+    SamplerPubLike,
+)
+from qiskit.primitives.containers.sampler_pub import SamplerPub
+from qiskit.primitives.primitive_job import PrimitiveJob
+...
+
+class CustomStatevectorSampler(BaseSamplerV2):
+
+    def __init__(self, *, default_shots: int = 1024, seed: np.random.Generator | int | None = None):
+        """
+        Args:
+            default_shots: The default shots for the sampler if not specified during run.
+            seed: The seed or Generator object for random number generation.
+                If None, a random seeded default RNG will be used.
+        """
+        self._default_shots = default_shots
+        self._seed = seed
+
+    ...
+
+    def run(
+        self, pubs: Iterable[SamplerPubLike], *, shots: int | None = None
+    ) -> PrimitiveJob[PrimitiveResult[SamplerPubResult]]:
+        if shots is None:
+            shots = self._default_shots
+        coerced_pubs = [SamplerPub.coerce(pub, shots) for pub in pubs]
+        if any(len(pub.circuit.cregs) == 0 for pub in coerced_pubs):
+            warnings.warn(
+                "One of your circuits has no output classical registers and so the result "
+                "will be empty. Did you mean to add measurement instructions?",
+                UserWarning,
+            )
+
+        job = PrimitiveJob(self._run, coerced_pubs)
+        job._submit()
+        return job
+
+    def _run(self, pubs: Iterable[SamplerPub]) -> PrimitiveResult[SamplerPubResult]:
+        results = [self._run_pub(pub) for pub in pubs]
+        return PrimitiveResult(results, metadata={"version": 2})
+
+    def _run_pub(self, pub: SamplerPub) -> SamplerPubResult:
+        circuit, qargs, meas_info = _preprocess_circuit(pub.circuit)
+        bound_circuits = pub.parameter_values.bind_all(circuit)
+        arrays = {
+            item.creg_name: np.zeros(
+                bound_circuits.shape + (pub.shots, item.num_bytes), dtype=np.uint8
+            )
+            for item in meas_info
+        }
+        for index, bound_circuit in enumerate(bound_circuits):
+
+            # ACCESS PROVIDER RESOURCE HERE
+            # in this case, we are showing an illustrative implementation
+            samples_array = ProviderResource.sample()
+
+            for item in meas_info:
+                ary = _samples_to_packed_array(samples_array, item.num_bits, item.qreg_indices)
+                arrays[item.creg_name][index] = ary
+
+        meas = {
+            item.creg_name: BitArray(arrays[item.creg_name], item.num_bits) for item in meas_info
+        }
+        return SamplerPubResult(
+            DataBin(**meas, shape=pub.shape),
+            metadata={"shots": pub.shots, "circuit_metadata": pub.circuit.metadata},
+        )    
+```
+
+The mechanics to implement a custom Estimator are analogous to the Sampler, but may require an extra pre/post-processing 
+step in the run method to extract expectation values from samples. In the case of the custom SatevectorEstimator that will be
+used as reference, we encapsulate thi computation in the `_statevector_from_circuit` utility function.
+
+
+``` python
+from .base import BaseEstimatorV2
+from .containers import DataBin, EstimatorPubLike, PrimitiveResult, PubResult
+from .containers.estimator_pub import EstimatorPub
+from .primitive_job import PrimitiveJob
+from .utils import _statevector_from_circuit
+...
+
+class CustomStatevectorEstimator(BaseEstimatorV2):
+
+    def __init__(
+        self, *, default_precision: float = 0.0, seed: np.random.Generator | int | None = None
+    ):
+        """
+        Args:
+            default_precision: The default precision for the estimator if not specified during run.
+            seed: The seed or Generator object for random number generation.
+                If None, a random seeded default RNG will be used.
+        """
+        self._default_precision = default_precision
+        self._seed = seed
+
+    ...
+
+    def run(
+        self, pubs: Iterable[EstimatorPubLike], *, precision: float | None = None
+    ) -> PrimitiveJob[PrimitiveResult[PubResult]]:
+        if precision is None:
+            precision = self._default_precision
+        coerced_pubs = [EstimatorPub.coerce(pub, precision) for pub in pubs]
+
+        job = PrimitiveJob(self._run, coerced_pubs)
+        job._submit()
+        return job
+
+    def _run(self, pubs: list[EstimatorPub]) -> PrimitiveResult[PubResult]:
+        return PrimitiveResult([self._run_pub(pub) for pub in pubs], metadata={"version": 2})
+
+    def _run_pub(self, pub: EstimatorPub) -> PubResult:
+        rng = np.random.default_rng(self._seed)
+        circuit = pub.circuit
+        observables = pub.observables
+        parameter_values = pub.parameter_values
+        precision = pub.precision
+        bound_circuits = parameter_values.bind_all(circuit)
+        bc_circuits, bc_obs = np.broadcast_arrays(bound_circuits, observables)
+        evs = np.zeros_like(bc_circuits, dtype=np.float64)
+        stds = np.zeros_like(bc_circuits, dtype=np.float64)
+        for index in np.ndindex(*bc_circuits.shape):
+            bound_circuit = bc_circuits[index]
+            observable = bc_obs[index]
+            final_state = _statevector_from_circuit(bound_circuit, rng)
+            paulis, coeffs = zip(*observable.items())
+            obs = SparsePauliOp(paulis, coeffs)
+            expectation_value = np.real_if_close(final_state.expectation_value(obs))
+            if precision != 0:
+                if not np.isreal(expectation_value):
+                    raise ValueError("Given operator is not Hermitian and noise cannot be added.")
+                expectation_value = rng.normal(expectation_value, precision)
+            evs[index] = expectation_value
+
+        data = DataBin(evs=evs, stds=stds, shape=evs.shape)
+        return PubResult(
+            data, metadata={"target_precision": precision, "circuit_metadata": pub.circuit.metadata}
+        )
+```
+