noisy_eval.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, transpile, execute\n",
    "from qiskit.providers.fake_provider import FakeCairo, FakeKolkata, FakeMontreal\n",
    "from qiskit_aer import AerSimulator\n",
    "from qiskit_aer.noise import NoiseModel\n",
    "from qiskit.result import LocalReadoutMitigator\n",
    "import pickle\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import numpy as np\n",
    "import qiskit.providers.aer.noise as noise\n",
    "import numpy as np\n",
    "import json\n",
    "from qiskit import pulse\n",
    "\n",
    "\n",
    "class exp_runner:\n",
    "    def __init__(self, model, amp_thresh=0.4):\n",
    "        self.model = (\n",
    "            model  # name of the error model, can be \"cairo\", \"kolkata\", \"montreal\"\n",
    "        )\n",
    "        self.amp_thresh = amp_thresh  # threshold for pauli string pruning\n",
    "        self.get_backend()\n",
    "        self.get_circ_and_obs()\n",
    "\n",
    "    def get_backend(self):\n",
    "        # get the backend and mitigator for the backend\n",
    "        if self.model == \"cairo\":\n",
    "            backend = FakeCairo()\n",
    "            with open(\"NoiseModel/fakecairo.pkl\", \"rb\") as f:\n",
    "                noise_model = pickle.load(f)\n",
    "            layout = [1,2,4,5,6,8,9,10,11,18,22,26]  # best layout that minimize noise\n",
    "        if self.model == \"kolkata\":\n",
    "            backend = FakeKolkata()\n",
    "            with open(\"NoiseModel/fakekolkata.pkl\", \"rb\") as f:\n",
    "                noise_model = pickle.load(f)\n",
    "            layout = [2, 3, 8, 9, 10, 12, 14, 16, 18, 19, 22, 26]\n",
    "        if self.model == \"montreal\":\n",
    "            backend = FakeMontreal()\n",
    "            with open(\"NoiseModel/fakemontreal.pkl\", \"rb\") as f:\n",
    "                noise_model = pickle.load(f)\n",
    "            layout = [0, 3, 5, 10, 11, 13, 14, 18, 21, 23, 24, 25]\n",
    "        noise_model = noise.NoiseModel.from_dict(noise_model)\n",
    "        simulator = AerSimulator(noise_model=noise_model)\n",
    "        mitigator = LocalReadoutMitigator(backend=backend, qubits=layout)\n",
    "        self.backend = backend\n",
    "        self.noise_model = noise_model\n",
    "        self.simulator = simulator\n",
    "        self.mitigator = mitigator\n",
    "        self.layout = layout\n",
    "        self.total_shot = 0  # a counter for the total number of shots used\n",
    "\n",
    "    def get_circ_and_obs(self):\n",
    "        # get the ansatz circuit and the shots for each pauli string\n",
    "        with open(\"Hamiltonian/OHhamiltonian_groups.json\", \"r\") as f:\n",
    "            groups = json.load(\n",
    "                f\n",
    "            )  # dict of dict, the first key is the superobservable, the second key is the sub-observable, the value is the coefficient\n",
    "        obs = groups.keys()\n",
    "        # a subobservable is a pauli string that can be generated by the original pauli string by replacing some of the pauli operators with I\n",
    "        subobs = groups\n",
    "        shots_dict = {}\n",
    "        # how many shots to use for each superobservable\n",
    "        total_amp = 0\n",
    "\n",
    "        for ob in obs:\n",
    "            amp = 0\n",
    "            for subob in subobs[ob]:\n",
    "                amp += subobs[ob][subob] ** 2\n",
    "            amp = amp**0.5\n",
    "            if amp > self.amp_thresh:\n",
    "                total_amp += amp\n",
    "        for ob in obs:\n",
    "            amp = 0\n",
    "            for subob in subobs[ob]:\n",
    "                amp += subobs[ob][subob] ** 2\n",
    "            amp = amp**0.5\n",
    "            if amp > self.amp_thresh:\n",
    "                shots_dict[ob] = np.floor(1800000 * amp / total_amp)\n",
    "            else:\n",
    "                shots_dict[ob] = 0\n",
    "        self.obs = obs\n",
    "        self.shots_dict = shots_dict\n",
    "        self.subobs = subobs\n",
    "        with open(\"ansatz.qasm\", \"r\") as f:\n",
    "            ansatz = f.read()\n",
    "            circ = QuantumCircuit.from_qasm_str(ansatz)\n",
    "        self.circ = circ\n",
    "\n",
    "    def diagonal_and_run(self, circ, obs, shots, seed):\n",
    "        # diagonalize the circuit and sample from it\n",
    "        new_circ = circ.copy()\n",
    "        for i in range(12):\n",
    "            if obs[i] == \"X\":\n",
    "                new_circ.h(i)\n",
    "            elif obs[i] == \"Y\":\n",
    "                new_circ.sdg(i)\n",
    "                new_circ.h(i)\n",
    "            elif obs[i] == \"Z\":\n",
    "                pass\n",
    "            else:\n",
    "                pass\n",
    "        new_circ.measure_all()\n",
    "        result = (\n",
    "            execute(\n",
    "                new_circ,\n",
    "                self.simulator,\n",
    "                initial_layout=self.layout,\n",
    "                shots=shots,\n",
    "                seed_simulator=seed,\n",
    "            )\n",
    "            .result()\n",
    "            .get_counts()\n",
    "        )\n",
    "        self.total_shot += shots\n",
    "        return result\n",
    "\n",
    "    def construct_diagonal_matrix(self, subob):\n",
    "        # construct the diagonal matrix for a sub-observable\n",
    "        diagonal = np.ones(1)\n",
    "        for i in range(12):\n",
    "            if subob[i] != \"I\":\n",
    "                diagonal = np.kron(np.array([1, -1]), diagonal)\n",
    "            else:\n",
    "                diagonal = np.kron(np.array([1, 1]), diagonal)\n",
    "        return diagonal\n",
    "\n",
    "    def calc_exp_val(self, circ, shots_adjustment=1, seed=0):\n",
    "        # calculate the expectation value of the ansatz circuit, ignoring identity pauli strings\n",
    "        exp_val = 0\n",
    "        for ob in self.obs:\n",
    "            shots = self.shots_dict[ob]\n",
    "            if shots == 0: # ignore pauli strings that have amplitude smaller than the threshold\n",
    "                continue\n",
    "            data = self.diagonal_and_run(circ, ob, shots * shots_adjustment, seed)\n",
    "            diagonal = np.zeros(2**12)\n",
    "            for subob in self.subobs[ob].keys():\n",
    "                diagonal += (\n",
    "                    self.construct_diagonal_matrix(subob) * self.subobs[ob][subob]\n",
    "                )\n",
    "            val = self.mitigator.expectation_value(data, diagonal=diagonal)[0]\n",
    "            exp_val += val\n",
    "        return exp_val\n",
    "\n",
    "    def run(self, seed):\n",
    "        # calculate the expectation value using ZNE, this is our method\n",
    "        self.total_shot = 0\n",
    "        circ = self.circ.copy()\n",
    "        circ = transpile(circ, self.simulator)\n",
    "\n",
    "        org_exp_val = self.calc_exp_val(circ, shots_adjustment=1.5 / 2, seed=seed) # use 75% of the shots for the original circuit\n",
    "\n",
    "        new_circ = self.circ.copy()\n",
    "        new_circ.barrier()\n",
    "        new_circ.extend(self.circ.inverse())\n",
    "        new_circ.barrier()\n",
    "        new_circ.extend(self.circ)\n",
    "        new_exp_val = self.calc_exp_val(new_circ, shots_adjustment=0.5 / 2, seed=seed) # use 25% of the shots for the new circuit\n",
    "        return (\n",
    "            1.5 * org_exp_val - 0.5 * new_exp_val - 51.1139071673799 + 4.36537496654537\n",
    "        )\n",
    "\n",
    "    def no_ZNE(self, seed):\n",
    "        # calculate the expectation value without ZNE, this is used for comparison\n",
    "        circ = self.circ.copy()\n",
    "        circ = transpile(circ, self.simulator)\n",
    "        org_exp_val = self.calc_exp_val(circ, shots_adjustment=1, seed=seed)\n",
    "        return org_exp_val - 51.1139071673799 + 4.36537496654537\n",
    "\n",
    "    def get_pulse_duration(self):\n",
    "        with pulse.build(self.backend) as my_program1:\n",
    "            with pulse.transpiler_settings(optimization_level=0):\n",
    "                pulse.call(self.circ)\n",
    "        print(f\"pulse duration:{my_program1.duration}\")\n",
    "\n",
    "    def get_total_shot(self):\n",
    "        print(f\"total shots:{self.total_shot}\")\n",
    "\n",
    "\n",
    "def calc_score(result):\n",
    "    print(f\"Score:{(1-np.abs(-74.38714627 - np.mean(result))/74.38714627)*100}\")\n",
    "\n",
    "\n",
    "def calc_variance(result):\n",
    "    print(f\"Variance:{np.std(result)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_list=[20, 21, 30, 33, 36, 42, 43, 55, 67, 170]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pulse duration:96\n",
      "total shots:1800000.0\n",
      "Score:99.89092044595634\n",
      "Variance:0.002954454427173438\n"
     ]
    }
   ],
   "source": [
    "model = \"cairo\"\n",
    "runner = exp_runner(model)\n",
    "result_list = []\n",
    "for seed in seed_list:\n",
    "    result_list.append(runner.run(seed))\n",
    "runner.get_pulse_duration()\n",
    "runner.get_total_shot()\n",
    "calc_score(result_list)\n",
    "calc_variance(result_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pulse duration:160\n",
      "total shots:1800000.0\n",
      "Score:99.88718060201029\n",
      "Variance:0.0021725095261954236\n"
     ]
    }
   ],
   "source": [
    "model = \"kolkata\"\n",
    "runner = exp_runner(model)\n",
    "result_list = []\n",
    "for seed in seed_list:\n",
    "    result_list.append(runner.run(seed))\n",
    "runner.get_pulse_duration()\n",
    "runner.get_total_shot()\n",
    "calc_score(result_list)\n",
    "calc_variance(result_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pulse duration:160\n",
      "total shots:1800000.0\n",
      "Score:99.88647132711903\n",
      "Variance:0.0027439212028582656\n"
     ]
    }
   ],
   "source": [
    "model = \"montreal\"\n",
    "runner = exp_runner(model)\n",
    "result_list = []\n",
    "for seed in seed_list:\n",
    "    result_list.append(runner.run(seed))\n",
    "runner.get_pulse_duration()\n",
    "runner.get_total_shot()\n",
    "calc_score(result_list)\n",
    "calc_variance(result_list)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "tq",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}