Merge pull request #617 from qiboteam/status

Implementing validation using Chi squared

doc/source/getting-started/index.rst

@@ -12,5 +12,6 @@ your quantum hardware.
     interface
     runcard
     autoruncard
+    validation
     protocols
     example

doc/source/getting-started/validation.rst

@@ -0,0 +1,23 @@
+How to add validation to your protocol?
+=======================================
+
+In Qibocal there is the possibility to add a validation step
+during the execution of your protocols.
+Here is an example of a runcard which validates the results through
+:math:`\chi^2`.
+
+.. code-block:: yaml
+
+   actions:
+    - id: t1
+      priority: 0
+      operation: t1
+      validator:
+        scheme: chi2
+        parameters:
+            chi2_max_value: 1
+      parameters:
+        ...
+
+The execution will be interrupted in this case if the :math:`\chi^2` exceeds
+`chi_max_value`.

src/qibocal/auto/execute.py

@@ -10,6 +10,7 @@
 from .graph import Graph
 from .history import History
 from .runcard import Id, Runcard
+from .status import Failure
 from .task import Qubits, Task
 
 
@@ -140,6 +141,10 @@ def run(self, mode):
                 mode=mode,
             )
             self.history.push(completed)
+            if isinstance(completed.status, Failure) and mode.name == "autocalibration":
+                log.warning("Stopping execution due to error in validation.")
+                yield task.uid
+                break
             self.head = self.next()
             update = self.update and task.update
             if (

src/qibocal/auto/runcard.py

@@ -10,6 +10,7 @@
 from qibolab.qubits import QubitId
 
 from .operation import OperationId
+from .validation import Validator
 
 Id = NewType("Id", str)
 """Action identifiers type."""
@@ -35,6 +36,8 @@ class Action:
     """Local qubits (optional)."""
     update: bool = True
     """Runcard update mechanism."""
+    validator: Optional[Validator] = None
+    """Define validation scheme and parameters."""
     parameters: Optional[dict[str, Any]] = None
     """Input parameters, either values or provider reference."""
 

src/qibocal/auto/status.py

@@ -38,5 +38,5 @@ class Normal(Status):
     """All green."""
 
 
-class Broken(Status):
+class Failure(Status):
     """Unrecoverable."""

src/qibocal/auto/task.py

@@ -26,7 +26,7 @@
     dummy_operation,
 )
 from .runcard import Action, Id
-from .status import Normal, Status
+from .status import Failure, Normal, Status
 
 MAX_PRIORITY = int(1e9)
 """A number bigger than whatever will be manually typed. But not so insanely big not to fit in a native integer."""
@@ -86,6 +86,20 @@ def operation(self):
 
         return Operation[self.action.operation].value
 
+    def validate(self, results: Results) -> Optional[Status]:
+        """Performs validation only if validator is provided."""
+
+        if self.action.validator is None:
+            return None
+        status = {}
+        for qubit in self.qubits:
+            status[qubit] = self.action.validator.__call__(results, qubit)
+        # exit if any of the qubit state is Failure
+        if any(isinstance(stat, Failure) for stat in status.values()):
+            return Failure()
+        else:
+            return Normal()
+
     @property
     def main(self):
         """Main node to be executed next."""
@@ -156,10 +170,9 @@ def run(
                 completed.data, completed.data_time = operation.acquisition(
                     parameters, platform=platform
                 )
-
         if mode.name in ["autocalibration", "fit"]:
             completed.results, completed.results_time = operation.fit(completed.data)
-
+            completed.status = self.validate(completed.results)
         return completed
 
 

src/qibocal/auto/validation.py

@@ -0,0 +1,30 @@
+"""Validation module."""
+from dataclasses import dataclass
+from typing import NewType, Optional, Union
+
+from qibolab.qubits import QubitId, QubitPairId
+
+from ..config import log
+from .operation import Results
+from .status import Status
+from .validators import VALIDATORS
+
+ValidatorId = NewType("ValidatorId", str)
+"""Identifier for validator object."""
+
+
+@dataclass
+class Validator:
+    """Generic validator object."""
+
+    scheme: ValidatorId
+    parameters: Optional[dict] = None
+
+    def __call__(
+        self, results: Results, qubit: Union[QubitId, QubitPairId, list[QubitId]]
+    ) -> Status:
+        log.info(
+            f"Performing validation in qubit {qubit} of {results.__class__.__name__} using {self.scheme} scheme."
+        )
+        validator = VALIDATORS[self.scheme]
+        return validator(results, qubit, **self.parameters)

src/qibocal/auto/validators/__init__.py

@@ -0,0 +1,4 @@
+from .chi2 import check_chi2
+
+VALIDATORS = {"chi2": check_chi2}
+"""Dict with available validators."""

src/qibocal/auto/validators/chi2.py

@@ -0,0 +1,37 @@
+"""Chi2 validation"""
+from typing import Union
+
+from qibolab.qubits import QubitId, QubitPairId
+
+from qibocal.config import raise_error
+
+from ..operation import Results
+from ..status import Failure, Normal
+
+CHI2_MAX = 0.05
+"""Max value for accepting fit result."""
+
+
+def check_chi2(
+    results: Results,
+    qubit: Union[QubitId, QubitPairId, list[QubitId]],
+    chi2_max_value=None,
+):
+    """Performs validation of results using chi2.
+
+    It assesses that chi2 is below the chi2_max_value threshold.
+
+    """
+
+    if chi2_max_value is None:
+        chi2_max_value = CHI2_MAX
+    try:
+        chi2 = getattr(results, "chi2")[qubit][0]
+        if chi2 < chi2_max_value:
+            return Normal()
+        else:
+            return Failure()
+    except AttributeError:
+        raise_error(
+            NotImplementedError, f"Chi2 validation not available for {type(results)}"
+        )

src/qibocal/protocols/characterization/coherence/spin_echo.py

@@ -1,4 +1,5 @@
-from dataclasses import dataclass
+from dataclasses import dataclass, field
+from typing import Optional
 
 import numpy as np
 import plotly.graph_objects as go
@@ -10,7 +11,7 @@
 from qibocal import update
 from qibocal.auto.operation import Parameters, Qubits, Results, Routine
 
-from ..utils import table_dict, table_html
+from ..utils import chi2_reduced, table_dict, table_html
 from . import t1
 from .utils import exp_decay, exponential_fit_probability
 
@@ -35,6 +36,10 @@ class SpinEchoResults(Results):
     """T2 echo for each qubit."""
     fitted_parameters: dict[QubitId, dict[str, float]]
     """Raw fitting output."""
+    chi2: Optional[dict[QubitId, tuple[float, Optional[float]]]] = field(
+        default_factory=dict
+    )
+    """Chi squared estimate mean value and error."""
 
 
 class SpinEchoData(t1.T1Data):
@@ -118,8 +123,19 @@ def _acquisition(
 def _fit(data: SpinEchoData) -> SpinEchoResults:
     """Post-processing for SpinEcho."""
     t2Echos, fitted_parameters = exponential_fit_probability(data)
+    chi2 = {
+        qubit: (
+            chi2_reduced(
+                data[qubit].prob,
+                exp_decay(data[qubit].wait, *fitted_parameters[qubit]),
+                data[qubit].error,
+            ),
+            np.sqrt(2 / len(data[qubit].prob)),
+        )
+        for qubit in data.qubits
+    }
 
-    return SpinEchoResults(t2Echos, fitted_parameters)
+    return SpinEchoResults(t2Echos, fitted_parameters, chi2)
 
 
 def _plot(data: SpinEchoData, qubit, fit: SpinEchoResults = None):
@@ -177,8 +193,8 @@ def _plot(data: SpinEchoData, qubit, fit: SpinEchoResults = None):
         fitting_report = table_html(
             table_dict(
                 qubit,
-                ["T2 Spin Echo [ns]"],
-                [fit.t2_spin_echo[qubit]],
+                ["T2 Spin Echo [ns]", "chi2 reduced"],
+                [fit.t2_spin_echo[qubit], fit.chi2[qubit]],
                 display_error=True,
             )
         )

src/qibocal/protocols/characterization/coherence/t1.py

@@ -1,4 +1,5 @@
 from dataclasses import dataclass, field
+from typing import Optional
 
 import numpy as np
 import numpy.typing as npt
@@ -12,7 +13,7 @@
 from qibocal import update
 from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
 
-from ..utils import table_dict, table_html
+from ..utils import chi2_reduced, table_dict, table_html
 from . import utils
 
 COLORBAND = "rgba(0,100,80,0.2)"
@@ -39,6 +40,10 @@ class T1Results(Results):
     """T1 for each qubit."""
     fitted_parameters: dict[QubitId, dict[str, float]]
     """Raw fitting output."""
+    chi2: Optional[dict[QubitId, tuple[float, Optional[float]]]] = field(
+        default_factory=dict
+    )
+    """Chi squared estimate mean value and error."""
 
 
 CoherenceProbType = np.dtype(
@@ -138,7 +143,19 @@ def _fit(data: T1Data) -> T1Results:
             y = p_0-p_1 e^{-x p_2}.
     """
     t1s, fitted_parameters = utils.exponential_fit_probability(data)
-    return T1Results(t1s, fitted_parameters)
+    chi2 = {
+        qubit: (
+            chi2_reduced(
+                data[qubit].prob,
+                utils.exp_decay(data[qubit].wait, *fitted_parameters[qubit]),
+                data[qubit].error,
+            ),
+            np.sqrt(2 / len(data[qubit].prob)),
+        )
+        for qubit in data.qubits
+    }
+
+    return T1Results(t1s, fitted_parameters, chi2)
 
 
 def _plot(data: T1Data, qubit, fit: T1Results = None):
@@ -191,12 +208,20 @@ def _plot(data: T1Data, qubit, fit: T1Results = None):
             )
         )
         fitting_report = table_html(
-            table_dict(qubit, ["T1 [ns]"], [fit.t1[qubit]], display_error=True)
+            table_dict(
+                qubit,
+                [
+                    "T1 [ns]",
+                    "chi2 reduced",
+                ],
+                [fit.t1[qubit], fit.chi2[qubit]],
+                display_error=True,
+            )
         )
     # last part
     fig.update_layout(
         showlegend=True,
-        xaxis_title="Time (ns)",
+        xaxis_title="Time [ns]",
         yaxis_title="Probability of State 1",
     )
 

src/qibocal/protocols/characterization/coherence/t2.py

@@ -1,4 +1,5 @@
-from dataclasses import dataclass
+from dataclasses import dataclass, field
+from typing import Optional
 
 import numpy as np
 import plotly.graph_objects as go
@@ -11,7 +12,7 @@
 from qibocal import update
 from qibocal.auto.operation import Parameters, Qubits, Results, Routine
 
-from ..utils import table_dict, table_html
+from ..utils import chi2_reduced, table_dict, table_html
 from . import t1, utils
 
 
@@ -35,6 +36,10 @@ class T2Results(Results):
     """T2 for each qubit (ns)."""
     fitted_parameters: dict[QubitId, dict[str, float]]
     """Raw fitting output."""
+    chi2: Optional[dict[QubitId, tuple[float, Optional[float]]]] = field(
+        default_factory=dict
+    )
+    """Chi squared estimate mean value and error."""
 
 
 class T2Data(t1.T1Data):
@@ -111,7 +116,18 @@ def _fit(data: T2Data) -> T2Results:
         y = p_0 - p_1 e^{-x p_2}.
     """
     t2s, fitted_parameters = utils.exponential_fit_probability(data)
-    return T2Results(t2s, fitted_parameters)
+    chi2 = {
+        qubit: (
+            chi2_reduced(
+                data[qubit].prob,
+                utils.exp_decay(data[qubit].wait, *fitted_parameters[qubit]),
+                data[qubit].error,
+            ),
+            np.sqrt(2 / len(data[qubit].prob)),
+        )
+        for qubit in data.qubits
+    }
+    return T2Results(t2s, fitted_parameters, chi2)
 
 
 def _plot(data: T2Data, qubit, fit: T2Results = None):
@@ -168,12 +184,19 @@ def _plot(data: T2Data, qubit, fit: T2Results = None):
             )
         )
         fitting_report = table_html(
-            table_dict(qubit, ["T2 [ns]"], [fit.t2[qubit]], display_error=True)
+            table_dict(
+                qubit,
+                [
+                    "T2 [ns]",
+                    "chi2 reduced",
+                ],
+                [fit.t2[qubit], fit.chi2[qubit]],
+                display_error=True,
+            )
         )
     fig.update_layout(
         showlegend=True,
-        uirevision="0",  # ``uirevision`` allows zooming while live plotting
-        xaxis_title="Time (ns)",
+        xaxis_title="Time [ns]",
         yaxis_title="Probability of State 1",
     )