{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "70e7454a",
   "metadata": {},
   "source": [
    "# Simulating Circuits Locally with the Neutral Atom Simulator\n",
    "\n",
    "The QDK includes a local **neutral atom device simulator** that models the behavior of neutral-atom quantum hardware, including qubit loss noise. This makes it useful for testing circuits and understanding noisy results before submitting to real hardware on Azure Quantum.\n",
    "\n",
    "This notebook shows how to:\n",
    "- Run circuits against the local simulator using both **Qiskit** (`NeutralAtomBackend`) and **Cirq** (`NeutralAtomSampler`)\n",
    "- Configure qubit loss noise via `NoiseConfig`\n",
    "- Interpret results that include loss markers, using the **accepted** (clean) and **raw** (all shots) result fields"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ac71a38",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "\n",
    "Install the `qdk` package with Qiskit and Cirq support:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc300bf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install \"qdk[qiskit,cirq]\" matplotlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "181fe979",
   "metadata": {},
   "source": [
    "After installing, restart the kernel if it was already running. Then verify imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12819c3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, transpile\n",
    "from qiskit.visualization import plot_histogram\n",
    "from qdk.qiskit import NeutralAtomBackend\n",
    "\n",
    "import cirq\n",
    "from qdk.cirq import NeutralAtomSampler\n",
    "\n",
    "from qdk.simulation import NoiseConfig\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import Counter\n",
    "\n",
    "print(\"Imports successful.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f173eb7c",
   "metadata": {},
   "source": [
    "## The circuit: 3-qubit GHZ state\n",
    "\n",
    "We use a 3-qubit GHZ circuit throughout this notebook. It produces the maximally entangled state:\n",
    "\n",
    "$$|\\text{GHZ}\\rangle = \\frac{1}{\\sqrt{2}}(|000\\rangle + |111\\rangle)$$\n",
    "\n",
    "This is a good noise demonstration circuit because:\n",
    "- Ideal results are simple: approximately 50% `000` and 50% `111`\n",
    "- Any other outcome or missing qubit clearly indicates noise\n",
    "- With 3 qubits, each shot has three independent chances to experience qubit loss\n",
    "\n",
    "We define both a Qiskit and a Cirq version of the same circuit below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad3f4e58",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_qubits = 3\n",
    "SHOTS = 500\n",
    "SEED = 42"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b7bbee4",
   "metadata": {},
   "source": [
    "---\n",
    "## Part A: Qiskit\n",
    "\n",
    "### Build and transpile the circuit\n",
    "\n",
    "The `NeutralAtomBackend` natively supports the gate set `{Rz, SX, CZ}`. Circuits written with higher-level gates (like `H` and `CX`) must first be transpiled into this native set.\n",
    "\n",
    "We transpile manually with `skip_transpilation=True` on the subsequent `run()` calls. This gives us full control over the decomposition and avoids the backend re-transpiling on each run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96082349",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build the high-level GHZ circuit\n",
    "ghz_qiskit = QuantumCircuit(n_qubits, n_qubits)\n",
    "ghz_qiskit.h(0)\n",
    "for i in range(n_qubits - 1):\n",
    "    ghz_qiskit.cx(i, i + 1)\n",
    "ghz_qiskit.measure(range(n_qubits), range(n_qubits))\n",
    "\n",
    "print(\"High-level circuit:\")\n",
    "print(ghz_qiskit.draw())\n",
    "\n",
    "# Transpile to the native gate set {Rz, SX, CZ}\n",
    "backend = NeutralAtomBackend()\n",
    "native_qiskit = transpile(ghz_qiskit, backend=backend)\n",
    "\n",
    "print(\"\\nNative gate circuit:\")\n",
    "print(native_qiskit.draw())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29def8fb",
   "metadata": {},
   "source": [
    "### Noiseless simulation\n",
    "\n",
    "Run the native circuit with no noise configured. We expect to see only the two ideal GHZ outcomes: `000` and `111`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c72e5935",
   "metadata": {},
   "outputs": [],
   "source": [
    "job_noiseless = backend.run(\n",
    "    native_qiskit,\n",
    "    shots=SHOTS,\n",
    "    seed=SEED,\n",
    "    skip_transpilation=True,\n",
    ")\n",
    "data_noiseless = job_noiseless.result().results[0].data\n",
    "\n",
    "print(\"Noiseless counts:\", dict(data_noiseless.counts))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "plot_histogram(dict(data_noiseless.counts), ax=ax)\n",
    "ax.set_title(\"Qiskit — Noiseless GHZ\")\n",
    "ax.yaxis.grid(True, linestyle=\"--\", alpha=0.7)\n",
    "ax.set_axisbelow(True)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc1ae90",
   "metadata": {},
   "source": [
    "### Configure qubit loss noise\n",
    "\n",
    "`NoiseConfig` controls the per-gate noise parameters. Here we set an 8% qubit-loss probability on `Rz` gates.\n",
    "\n",
    "Since the `H` gate decomposes to `Rz` gates in the native set, every qubit passes through at least one `Rz` and has a chance of being lost before measurement. When a qubit is lost, the simulator records `\"-\"` in its bitstring position rather than `\"0\"` or `\"1\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70749b30",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = NoiseConfig()\n",
    "noise.rz.loss = 0.08  # 8% qubit-loss probability per Rz gate"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae8361a6",
   "metadata": {},
   "source": [
    "### Noisy simulation and reading results\n",
    "\n",
    "The result data (`job.result().results[0].data`) exposes two parallel sets of fields:\n",
    "\n",
    "| Field | What it contains |\n",
    "|---|---|\n",
    "| **`counts`** | Bitstring → shot count, accepted shots only (no `\"-\"`) |\n",
    "| **`probabilities`** | Bitstring → empirical probability, accepted shots only |\n",
    "| **`memory`** | Per-shot bitstring list, accepted shots only |\n",
    "| **`raw_counts`** | Bitstring → shot count, all shots (loss shots have `\"-\"`) |\n",
    "| **`raw_probabilities`** | Bitstring → empirical probability, all shots |\n",
    "| **`raw_memory`** | Per-shot bitstring list, all shots |\n",
    "\n",
    "Use the **accepted** fields for downstream analysis. Use the **raw** fields to inspect or quantify loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bb7d3f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "job_noisy = backend.run(\n",
    "    native_qiskit,\n",
    "    shots=SHOTS,\n",
    "    noise=noise,\n",
    "    seed=SEED,\n",
    "    skip_transpilation=True,\n",
    ")\n",
    "data_noisy = job_noisy.result().results[0].data\n",
    "\n",
    "accepted = dict(data_noisy.counts)\n",
    "raw      = dict(data_noisy.raw_counts)\n",
    "\n",
    "total_raw      = sum(raw.values())\n",
    "total_accepted = sum(accepted.values())\n",
    "total_lost     = total_raw - total_accepted\n",
    "\n",
    "print(f\"Total shots : {total_raw}\")\n",
    "print(f\"Accepted    : {total_accepted}  ({100 * total_accepted / total_raw:.1f}%)\")\n",
    "print(f\"Lost        : {total_lost}  ({100 * total_lost / total_raw:.1f}%)\")\n",
    "print()\n",
    "print(\"Accepted counts:\", accepted)\n",
    "print(\"Raw counts     :\", raw)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63b0aa3e",
   "metadata": {},
   "source": [
    "### Visualize: accepted vs raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9cead95",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n",
    "\n",
    "plot_histogram(accepted, ax=axes[0])\n",
    "plot_histogram(raw,      ax=axes[1])\n",
    "\n",
    "for ax, title in zip(axes, [\n",
    "    \"Qiskit — Accepted shots (loss filtered out)\",\n",
    "    \"Qiskit — Raw shots (loss bitstrings included)\",\n",
    "]):\n",
    "    ax.set_title(title)\n",
    "    ax.yaxis.grid(True, linestyle=\"--\", alpha=0.7)\n",
    "    ax.set_axisbelow(True)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "700f305b",
   "metadata": {},
   "source": [
    "---\n",
    "## Part B: Cirq\n",
    "\n",
    "### Build the circuit\n",
    "\n",
    "The `NeutralAtomSampler` accepts standard Cirq circuits directly. It internally converts them to OpenQASM 3.0 before simulating, so no manual transpilation step is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3448c0a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "qubits = cirq.LineQubit.range(n_qubits)\n",
    "\n",
    "ghz_cirq = cirq.Circuit(\n",
    "    cirq.H(qubits[0]),\n",
    "    *[cirq.CNOT(qubits[i], qubits[i + 1]) for i in range(n_qubits - 1)],\n",
    "    cirq.measure(*qubits, key=\"result\"),\n",
    ")\n",
    "\n",
    "print(ghz_cirq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3ec4b8a",
   "metadata": {},
   "source": [
    "### Noiseless simulation\n",
    "\n",
    "Run with no noise and confirm only `000` and `111` appear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b2d36ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "sampler = NeutralAtomSampler()\n",
    "cirq_result_noiseless = sampler.run(ghz_cirq, repetitions=SHOTS)\n",
    "\n",
    "# histogram() returns a Counter of integer bitmask → count; format as binary strings\n",
    "cirq_noiseless_counts = {\n",
    "    format(k, f\"0{n_qubits}b\"): v\n",
    "    for k, v in cirq_result_noiseless.histogram(key=\"result\").items()\n",
    "}\n",
    "print(\"Noiseless counts:\", cirq_noiseless_counts)\n",
    "\n",
    "def plot_bar(ax, counts, title):\n",
    "    labels = sorted(counts.keys())\n",
    "    values = [counts[k] for k in labels]\n",
    "    bars = ax.bar(labels, values)\n",
    "    ax.set_title(title)\n",
    "    ax.set_xlabel(\"Outcome\")\n",
    "    ax.set_ylabel(\"Count\")\n",
    "    ax.tick_params(axis=\"x\", rotation=45)\n",
    "    ax.yaxis.grid(True, linestyle=\"--\", alpha=0.7)\n",
    "    ax.set_axisbelow(True)\n",
    "    for bar, value in zip(bars, values):\n",
    "        ax.text(\n",
    "            bar.get_x() + bar.get_width() / 2,\n",
    "            bar.get_height(), str(value),\n",
    "            ha=\"center\", va=\"bottom\", fontsize=9,\n",
    "        )\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "plot_bar(ax, cirq_noiseless_counts, \"Cirq — Noiseless GHZ\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2318fdd",
   "metadata": {},
   "source": [
    "### Noisy simulation and reading results\n",
    "\n",
    "Pass the same `NoiseConfig` to `NeutralAtomSampler`. The result object exposes:\n",
    "\n",
    "| Field | What it contains |\n",
    "|---|---|\n",
    "| **`result.measurements[key]`** | NumPy int8 array of accepted shots only (no `\"-\"`), shape `(accepted_shots, n_qubits)` |\n",
    "| **`result.raw_measurements()[key]`** | String array of all shots, `\"-\"` where a qubit was lost, same shape |\n",
    "| **`result.raw_shots`** | The original shot objects as returned by the simulator |\n",
    "\n",
    "Use `result.measurements` for analysis. Use `result.raw_measurements()` and `result.raw_shots` to inspect loss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f067bfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "noisy_sampler = NeutralAtomSampler(noise=noise, seed=SEED)\n",
    "cirq_result_noisy = noisy_sampler.run(ghz_cirq, repetitions=SHOTS)\n",
    "\n",
    "# Accepted: int8 array, shape = (accepted_shots, n_qubits)\n",
    "accepted_arr = cirq_result_noisy.measurements[\"result\"]\n",
    "\n",
    "# Raw: string array, shape = (total_shots, n_qubits), \"-\" for lost qubits\n",
    "raw_arr = cirq_result_noisy.raw_measurements()[\"result\"]\n",
    "\n",
    "cirq_total    = raw_arr.shape[0]\n",
    "cirq_accepted = accepted_arr.shape[0]\n",
    "cirq_lost     = cirq_total - cirq_accepted\n",
    "\n",
    "print(f\"Total shots : {cirq_total}\")\n",
    "print(f\"Accepted    : {cirq_accepted}  ({100 * cirq_accepted / cirq_total:.1f}%)\")\n",
    "print(f\"Lost        : {cirq_lost}  ({100 * cirq_lost / cirq_total:.1f}%)\")\n",
    "\n",
    "def arr_to_counts(arr):\n",
    "    return Counter(\"\".join(row) for row in arr)\n",
    "\n",
    "cirq_accepted_counts = arr_to_counts(accepted_arr.astype(str))\n",
    "cirq_raw_counts      = arr_to_counts(raw_arr)\n",
    "\n",
    "print(\"\\nAccepted counts:\", dict(cirq_accepted_counts))\n",
    "print(\"Raw counts     :\", dict(cirq_raw_counts))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "486fa6d4",
   "metadata": {},
   "source": [
    "### Visualize: accepted vs raw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4898a60",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n",
    "plot_bar(axes[0], dict(cirq_accepted_counts), \"Cirq — Accepted shots (loss filtered out)\")\n",
    "plot_bar(axes[1], dict(cirq_raw_counts),      \"Cirq — Raw shots (loss bitstrings included)\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06b1d027",
   "metadata": {},
   "source": [
    "## Notes\n",
    "\n",
    "- **`NoiseConfig` gate coverage**: Qubit-loss probability can be set independently per gate type (e.g. `noise.rz.loss`, `noise.sx.loss`, `noise.cz.loss`). You can mix different rates to model realistic hardware calibration profiles.\n",
    "- **Why `skip_transpilation=True` (Qiskit only)**: Passing `skip_transpilation=True` tells the backend that the circuit is already in the native gate set. If you omit it, the backend will transpile automatically, but transpiling once and reusing saves time when running many shots or noise configurations.\n",
    "- **Cirq result arrays**: Cirq's `measurements` and `raw_measurements()` return 2D NumPy arrays with one row per shot and one column per qubit. The `arr_to_counts` helper above joins each row into a single bitstring (e.g. `[\"1\", \"1\", \"1\"]` → `\"111\"`) to make the distribution easy to inspect and visualize.\n",
    "- **Next steps**: Once your circuit behaves as expected locally, submit it to real neutral-atom hardware on Azure Quantum using `cirq_submission_to_azure.ipynb` or `qiskit_submission_to_azure.ipynb`.\n"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
