notebooks/stage-two-analysis.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6bed6628",
   "metadata": {},
   "source": [
    "# TGP Stage 2 analysis\n",
    "\n",
    "This notebook performs an interactive analysis of one of the **Stage 2** datasets.\n",
    "It starts by loading the raw data and proceeds with pre-processing and analysis steps, where each step is visualized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7f7bb55",
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import annotations\n",
    "\n",
    "from pathlib import Path\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import xarray as xr\n",
    "import tgp\n",
    "\n",
    "tgp.plot.set_mpl_rc()\n",
    "\n",
    "folder = Path(\"..\") / \"data\"\n",
    "\n",
    "print(\"Using tgp code version:\", tgp.__version__)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a4f7fa7b",
   "metadata": {},
   "source": [
    "# Analysis on experimental data\n",
    "\n",
    "First, we select the dataset we want to analyze. The dictionary below contains all of the datasets that are provided.\n",
    "\n",
    "Change the `selected = \"A1\"` parameter to e.g., `selected = \"B\"` to analyze the data of device B instead of device A1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9f77971",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import NamedTuple\n",
    "\n",
    "\n",
    "class Parameters(NamedTuple):\n",
    "    lockin_1: str\n",
    "    lockin_2: str\n",
    "    fridge: str\n",
    "    drop_indices: list[int] | tuple[list[int], list[int]]\n",
    "    max_bias_index: int | tuple[int, int]\n",
    "    selected_cutter: int\n",
    "    selected_clusters: list[int]\n",
    "    gap_threshold_high: float = 70e-3\n",
    "    phase_shift_left: float = 0.0\n",
    "    phase_shift_right: float = 0.0\n",
    "\n",
    "\n",
    "parameters_dict = {\n",
    "    \"A1\": Parameters(  # ✅ passed TGP\n",
    "        lockin_1=\"lockin_1\",\n",
    "        lockin_2=\"lockin_2\",\n",
    "        fridge=\"deviceA1.yaml\",\n",
    "        drop_indices=[0, 1, 2],\n",
    "        max_bias_index=-2,\n",
    "        selected_cutter=0,\n",
    "        selected_clusters=[1],\n",
    "        phase_shift_left=-3.3,\n",
    "        phase_shift_right=-5.57,\n",
    "    ),\n",
    "    \"A2\": Parameters(  # ✅ passed TGP\n",
    "        lockin_1=\"mfli_5510\",\n",
    "        lockin_2=\"mfli_5602\",\n",
    "        fridge=\"deviceA2.yaml\",\n",
    "        drop_indices=[0, 1, 2],\n",
    "        max_bias_index=-2,\n",
    "        selected_cutter=1,\n",
    "        selected_clusters=[1],\n",
    "    ),\n",
    "    \"A3\": Parameters(  # ✅ passed TGP\n",
    "        lockin_1=\"mfli_5510\",\n",
    "        lockin_2=\"mfli_5602\",\n",
    "        fridge=\"deviceA2.yaml\",\n",
    "        drop_indices=[],\n",
    "        max_bias_index=-1,\n",
    "        selected_cutter=1,\n",
    "        selected_clusters=[1],\n",
    "    ),\n",
    "    \"B\": Parameters(  # ✅ passed TGP\n",
    "        lockin_1=\"lockin_1\",\n",
    "        lockin_2=\"lockin_2\",\n",
    "        fridge=\"deviceB.yaml\",\n",
    "        drop_indices=[0, 1, 2, 3],\n",
    "        max_bias_index=-4,\n",
    "        selected_cutter=1,\n",
    "        selected_clusters=[1],\n",
    "        gap_threshold_high=100e-3,\n",
    "        phase_shift_left=-3.78,\n",
    "        phase_shift_right=-6.89,\n",
    "    ),\n",
    "    \"C\": Parameters(  # ✅ passed TGP\n",
    "        lockin_1=\"mfli_5583\",\n",
    "        lockin_2=\"mfli_5591\",\n",
    "        fridge=\"deviceC.yaml\",\n",
    "        drop_indices=([0, 1, 2, 3, 4, 5], []),\n",
    "        max_bias_index=(-4, -1),\n",
    "        selected_cutter=0,\n",
    "        selected_clusters=[1],\n",
    "    ),\n",
    "    \"D\": Parameters(  # ✅ passed TGP\n",
    "        lockin_1=\"mfli_4909\",\n",
    "        lockin_2=\"mfli_5654\",\n",
    "        fridge=\"deviceD.yaml\",\n",
    "        drop_indices=[0, 1],\n",
    "        max_bias_index=-4,\n",
    "        selected_cutter=0,\n",
    "        selected_clusters=[1, 2],\n",
    "    ),\n",
    "    \"E\": Parameters(  # ❌ did not pass TGP\n",
    "        lockin_1=\"mfli_5591\",\n",
    "        lockin_2=\"mfli_5583\",\n",
    "        fridge=\"deviceE.yaml\",\n",
    "        drop_indices=[0],\n",
    "        max_bias_index=-1,\n",
    "        selected_cutter=1,\n",
    "        selected_clusters=[],\n",
    "    ),\n",
    "}\n",
    "\n",
    "# Select the dataset we want\n",
    "selected = \"A1\"\n",
    "p = parameters_dict[selected]\n",
    "fname_left = f\"device{selected}_stage2_left.nc\"\n",
    "fname_right = f\"device{selected}_stage2_right.nc\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "15044174",
   "metadata": {},
   "source": [
    "We load the raw data and prepare it for the `tgp` code's analysis.\n",
    "We call the `prepare` function on the loaded raw datasets to validate that it contains the correct dimensions and variables, it also performs renames (if necessary) to account for different naming schemes in different fridges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f1e0e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_left = xr.load_dataset(folder / \"experimental\" / fname_left, engine=\"h5netcdf\")\n",
    "ds_right = xr.load_dataset(folder / \"experimental\" / fname_right, engine=\"h5netcdf\")\n",
    "ds_left = tgp.prepare.prepare(ds_left)\n",
    "ds_right = tgp.prepare.prepare(ds_right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7258db9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "if selected == \"D\":  # Device D is for a single cutter pair only\n",
    "    # the code expects this dimension so we add it.\n",
    "    ds_left = ds_left.expand_dims(cutter_pair_index=[0])\n",
    "    ds_right = ds_right.expand_dims(cutter_pair_index=[0])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "a4215704",
   "metadata": {},
   "source": [
    "We can view the renamed raw data with the plot below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c406493",
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the data (before corrections)\n",
    "\n",
    "print(f\"Choose V={ds_left.V.values}\")\n",
    "print(f\"Choose cutter_pair_index={ds_left.cutter_pair_index.values}\")\n",
    "\n",
    "tgp.plot.two.plot_data(\n",
    "    ds_left,\n",
    "    ds_right,\n",
    "    V=ds_left.V.median(),\n",
    "    cutter_pair_index=p.selected_cutter,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cb8accc2",
   "metadata": {},
   "source": [
    "We apply corrections for the circuit effects in three-terminal electrical transport measurements arising from finite line impedances (as reported in [arXiv:2104.02671](https://arxiv.org/abs/2104.02671)).\n",
    "The `correct_frequencies` functions corrects the spurious voltage divider effects.\n",
    "Additionally, the `correct_bias` function ensures that zero-bias is at the correct point.\n",
    "This is required because zero-bias in the experiment might be slightly offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e01a304",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply corrections\n",
    "from tgp.frequency_correction import correct_frequencies\n",
    "\n",
    "ds_left, ds_right = correct_frequencies(\n",
    "    ds_left,\n",
    "    ds_right,\n",
    "    lockin_left=p.lockin_1,  # key in the snapshot\n",
    "    lockin_right=p.lockin_2,\n",
    "    fridge_parameters=folder / \"fridge\" / p.fridge,\n",
    "    phase_shift_left=p.phase_shift_left,\n",
    "    phase_shift_right=p.phase_shift_right,\n",
    "\n",
    ")\n",
    "\n",
    "ds_left, ds_right = tgp.prepare.correct_bias(\n",
    "    ds_left,\n",
    "    ds_right,\n",
    "    drop_indices=p.drop_indices,\n",
    "    max_bias_index=p.max_bias_index,\n",
    "    norm=1e3,\n",
    "    method=\"manual\",\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "64bdec22",
   "metadata": {},
   "source": [
    "We can view the conductance data again with the plot we used above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28a8ed01",
   "metadata": {},
   "outputs": [],
   "source": [
    "# View the data (after corrections)\n",
    "\n",
    "print(f\"Choose V={ds_left.V.values}\")\n",
    "print(f\"Choose cutter_pair_index={ds_left.cutter_pair_index.values}\")\n",
    "\n",
    "tgp.plot.two.plot_data(\n",
    "    ds_left,\n",
    "    ds_right,\n",
    "    V=ds_left.V.median(),  # choose the voltage here\n",
    "    cutter_pair_index=p.selected_cutter,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "3b9d66e3",
   "metadata": {},
   "source": [
    "We can now perform the first analysis step, which is to extract the value of the gap.\n",
    "This algorithm is based on thresholding the filtered anti-symmetrized part of the nonlocal conductance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9ad6c8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform the gap extration, optionally tweak the parameters by passing them\n",
    "ds_left, ds_right = tgp.two.extract_gap(ds_left, ds_right)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7e012152",
   "metadata": {},
   "source": [
    "We can see the results of the gap extraction using the plot below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5a285d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results of the gap extration\n",
    "\n",
    "print(f\"Choose V={ds_left.V.values}\")\n",
    "print(f\"Choose cutter_pair_index={ds_left.cutter_pair_index.values}\")\n",
    "\n",
    "V = ds_left.V.median()\n",
    "\n",
    "tgp.plot.two.plot_gap_extraction(ds_left, ds_right, p.selected_cutter, V)\n",
    "tgp.plot.two.plot_extracted_gap(ds_left, ds_right, p.selected_cutter)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "110dc893",
   "metadata": {},
   "source": [
    "In addition to learning in which region the phase diagram is gapped, we need to know where in `(B, V)` space zero-bias peaks (ZPB)s occur.\n",
    "The `zbp_dataset_derivative` function returns a dataset with a combined left and right gap and ZBPs, without the `left_bias` or `right_bias` dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "027514cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "zbp_ds = tgp.two.zbp_dataset_derivative(\n",
    "    ds_left, ds_right, zbp_probability_threshold=0.6\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "173629ea",
   "metadata": {},
   "source": [
    "To go from a gap array with numerical values to a boolean array where in `(B, V)` space the spectrum is gapped, we apply an upper and lower threshold."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "91ddd2a2",
   "metadata": {},
   "source": [
    "We can set them with the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cd72895",
   "metadata": {},
   "outputs": [],
   "source": [
    "tgp.two.set_gap_threshold(\n",
    "    zbp_ds,\n",
    "    threshold_low=10e-3,\n",
    "    threshold_high=p.gap_threshold_high,\n",
    ")\n",
    "\n",
    "tgp.plot.zbp.plot_gapped_zbp(zbp_ds, p.selected_cutter)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "745b31ae",
   "metadata": {},
   "source": [
    "Finally, we cluster the gapped ZBP array and score the clusters with the `cluster_and_score` function.\n",
    "This action is performed on every `cutter_pair_value` slice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42b4c9f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform analysis on all cutter_pair_index values\n",
    "zbp_ds_done = tgp.two.cluster_and_score(zbp_ds, min_cluster_size=7)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9d3dae19",
   "metadata": {},
   "source": [
    "We can plot the results of the analysis with the plot below.\n",
    "\n",
    "We can see the probabilities of the ZPBs on the left and right, and their joint probabilities (top row).\n",
    "Below we see the gap (d), the gapped ZPBs (e), and the clusters of the gapped ZBPs (f).\n",
    "Then in the lowest row, we plot the region of interest (ROI) and the gap inside that region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f879fa5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Choose one of cutter_pair_index={zbp_ds_done.cutter_pair_index.values}\")\n",
    "print(f\"Choose one of zbp_cluster_number={zbp_ds_done.zbp_cluster_number.values}\")\n",
    "\n",
    "if \"cutter_pair_index\" in zbp_ds_done.dims:\n",
    "    sel = zbp_ds_done.sel(cutter_pair_index=p.selected_cutter)\n",
    "else:\n",
    "    sel = zbp_ds_done\n",
    "tgp.plot.zbp.plot_probability_and_clusters(sel)\n",
    "if p.selected_clusters:\n",
    "    tgp.plot.zbp.plot_region_of_interest_2(sel, p.selected_clusters[0])\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d98f6da0",
   "metadata": {},
   "source": [
    "## Manually select cutter value\n",
    "\n",
    "The `cluster_and_score` function used above performs the entire analysis on a per `cutter_pair_index` basis and joins the results into an array with the `cutter_pair_index` dimension.\n",
    "Alternatively, we can select a single cutter value and use the analysis and plotting functions directly.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "The functions below assume a single cutter value.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b854784f",
   "metadata": {},
   "outputs": [],
   "source": [
    "zbp_ds_sel = (\n",
    "    zbp_ds.isel(cutter_pair_index=p.selected_cutter)\n",
    "    if \"cutter_pair_index\" in zbp_ds.dims\n",
    "    else zbp_ds\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "114d7d74",
   "metadata": {},
   "outputs": [],
   "source": [
    "tgp.plot.zbp.plot_left_right_zbp_probability(zbp_ds)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f2a3fc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "tgp.plot.zbp.plot_joined_probability(zbp_ds)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21851b8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tgp.plot.zbp.plot_clusters_zbp(zbp_ds_sel)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84a063de",
   "metadata": {},
   "outputs": [],
   "source": [
    "tgp.plot.zbp.plot_probability_and_clusters(zbp_ds_sel)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce846380",
   "metadata": {},
   "outputs": [],
   "source": [
    "if p.selected_clusters:\n",
    "    tgp.plot.zbp.plot_region_of_interest_2(zbp_ds_sel, zbp_cluster_number=p.selected_clusters[0])\n",
    "    plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e87f4951",
   "metadata": {},
   "source": [
    "## Paper plots of the selected cluster\n",
    "\n",
    "By selecting a single `cutter_pair_index` we can produce the same plot as in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a2990af",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tgp.plot.paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aea97a03",
   "metadata": {},
   "outputs": [],
   "source": [
    "tgp.plot.paper.plot_stage2_diagram(\n",
    "    zbp_ds,\n",
    "    cutter_value=p.selected_cutter,\n",
    "    zbp_cluster_numbers=p.selected_clusters,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "664c9c3b",
   "metadata": {},
   "source": [
    "One has to select same `selected_plunger` parameter as in the paper figures notebook to reproduce the same plot. We choose the plunger to be in the center of the selected cluster here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "454a52fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "if \"cutter_pair_index\" not in zbp_ds.dims:\n",
    "    zbp_ds = zbp_ds.expand_dims(cutter_pair_index=[0])\n",
    "\n",
    "if p.selected_clusters:\n",
    "    V = zbp_ds.sel(\n",
    "        zbp_cluster_number=p.selected_clusters[0], cutter_pair_index=p.selected_cutter\n",
    "    ).cluster_V_center.item()\n",
    "else:\n",
    "    raise RuntimeError(\"Please manually set the V value.\")\n",
    "\n",
    "tgp.plot.paper.plot_conductance_waterfall(\n",
    "    zbp_ds_sel,\n",
    "    ds_left,\n",
    "    ds_right,\n",
    "    selected_cutter=p.selected_cutter,\n",
    "    selected_plunger=V,\n",
    "    zbp_cluster_numbers=p.selected_clusters,\n",
    "    bias_max=0.06,\n",
    "    labels=\"abcd\",\n",
    ")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}