{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Quara's standard tomography features from Forest (2-qubit)\n",
    "\n",
    "Quara supports object conversions for executing standard quantum tomography from Forest SDK. Here we briefly explain how to perform quantum tomography by using Quara and measurements obtained from quantum virtual machine of Forest SDK.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "\n",
    "# quara\n",
    "from quara.interface.forest.api import (\n",
    "    generate_pauli_strings_from_povm_name,\n",
    "    calc_empi_dist_from_observables\n",
    ")\n",
    "from quara.objects.composite_system_typical import generate_composite_system\n",
    "from quara.objects.povm_typical import generate_povm_from_name\n",
    "from quara.protocol.qtomography.standard.standard_qst import StandardQst\n",
    "from quara.protocol.qtomography.standard.linear_estimator import LinearEstimator\n",
    "\n",
    "# pyquil\n",
    "from pyquil import get_qc, Program\n",
    "from pyquil.gates import H, CNOT\n",
    "from pyquil.paulis import PauliTerm\n",
    "from pyquil.experiment import (\n",
    "    Experiment,\n",
    "    ExperimentSetting,\n",
    "    zeros_state,\n",
    ")\n",
    "from pyquil.operator_estimation import measure_observables\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we define a program for state preparation for tomographic experiment. Then, we create a function `obtain_expectations_for_qst()` which is a function that performs measuments based on given pauli operators. Pauli operators are operators that can be performed as measurements in Forest SDK. More details will be explained in the next cell. Make sure you have started `quilc` and `qvm` on your other consoles. See [this page](https://pyquil-docs.rigetti.com/en/stable/start.html#setting-up-server-mode-for-pyquil) for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example for 2 qubits system\n",
    "qc = get_qc(\"2q-qvm\")\n",
    "qubits = [0, 1]\n",
    "\n",
    "# define initialization of the quantum system\n",
    "# creating a state called \"bell phi plus state\"\n",
    "num_shots = 10000\n",
    "p = Program()\n",
    "p += H(qubits[0])\n",
    "p += CNOT(qubits[0], qubits[1])\n",
    "p.wrap_in_numshots_loop(num_shots)\n",
    "\n",
    "def obtain_expectations_for_qst(qc, program, pauli_strings):\n",
    "    settings = []\n",
    "    for pauli_str in pauli_strings:\n",
    "        out_operator = PauliTerm.from_list(list(zip(pauli_str, qubits)))\n",
    "        settings.append(ExperimentSetting(zeros_state(qubits), out_operator))\n",
    "    tomo_experiment = Experiment(settings, program)\n",
    "    expectations = []\n",
    "    for pauli_str, res in zip(\n",
    "        pauli_strings,\n",
    "        measure_observables(\n",
    "            qc,\n",
    "            tomo_experiment,\n",
    "        ),\n",
    "    ):\n",
    "        if res.raw_expectation is None:\n",
    "            # This is the result for II...I operator\n",
    "            expectations.append(1.0)\n",
    "        else:\n",
    "            expectations.append(res.raw_expectation)\n",
    "        print(f\"Raw expectation of {pauli_str}: {expectations[-1]}\")\n",
    "    return expectations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define a composite system and POVM objects which will make a tomographically complete measuring experiment. Then we calculate a set of desired pauli mesruments by using a pre-difned function named `generate_pauli_strings_from_povm_name()` in Quara for each POVM. In this example, we create a dictionary for POVM objects and Pauli operator names which both are indexed with `povm_name`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample of Quara tester object for povm_name: x_x\n",
      "Type:\n",
      "Povm\n",
      "\n",
      "Dim:\n",
      "4\n",
      "\n",
      "Number of outcomes:\n",
      "4\n",
      "\n",
      "Vecs:\n",
      "[[ 0.5  0.5  0.   0.   0.5  0.5  0.   0.   0.   0.   0.   0.   0.   0.\n",
      "   0.   0. ]\n",
      " [ 0.5 -0.5  0.   0.   0.5 -0.5  0.   0.   0.   0.   0.   0.   0.   0.\n",
      "   0.   0. ]\n",
      " [ 0.5  0.5  0.   0.  -0.5 -0.5  0.   0.   0.   0.   0.   0.   0.   0.\n",
      "   0.   0. ]\n",
      " [ 0.5 -0.5  0.   0.  -0.5  0.5  0.   0.   0.   0.   0.   0.   0.   0.\n",
      "   0.   0. ]]\n",
      "\n",
      "Pauli strings that will be used in this tomographical experiment:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'x_x': ['XX', 'XI', 'IX', 'II'],\n",
       " 'x_y': ['XY', 'XI', 'IY', 'II'],\n",
       " 'x_z': ['XZ', 'XI', 'IZ', 'II'],\n",
       " 'y_x': ['YX', 'YI', 'IX', 'II'],\n",
       " 'y_y': ['YY', 'YI', 'IY', 'II'],\n",
       " 'y_z': ['YZ', 'YI', 'IZ', 'II'],\n",
       " 'z_x': ['ZX', 'ZI', 'IX', 'II'],\n",
       " 'z_y': ['ZY', 'ZI', 'IY', 'II'],\n",
       " 'z_z': ['ZZ', 'ZI', 'IZ', 'II']}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_sys = generate_composite_system(\"qubit\", 2)\n",
    "povms = {}\n",
    "pauli_strings = {}\n",
    "povm_names = [\"_\".join(i) for i in product(\"xyz\", repeat=2)]\n",
    "for povm_name in povm_names:\n",
    "    povms[povm_name] = generate_povm_from_name(povm_name, c_sys)\n",
    "    pauli_strings[povm_name] = generate_pauli_strings_from_povm_name(povm_name)\n",
    "print(f\"Sample of Quara tester object for povm_name: {povm_names[0]}\")\n",
    "print(povms[povm_names[0]])\n",
    "print(\"\\nPauli strings that will be used in this tomographical experiment:\")\n",
    "pauli_strings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we perform projective measurements based on Pauli operators and `obtain_expectations_for_qst()` which we defined in previous cells.\n",
    "Observed expectations are also stored as dictionary and indexed with `povm_name`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Raw expectation of XX: 1.0\n",
      "Raw expectation of XI: -0.0122\n",
      "Raw expectation of IX: -0.0052\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of XY: 0.0194\n",
      "Raw expectation of XI: 0.0054\n",
      "Raw expectation of IY: 0.0002\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of XZ: 0.0218\n",
      "Raw expectation of XI: -0.0052\n",
      "Raw expectation of IZ: -0.0022\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of YX: 0.0092\n",
      "Raw expectation of YI: -0.016\n",
      "Raw expectation of IX: -0.0028\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of YY: -1.0\n",
      "Raw expectation of YI: 0.0172\n",
      "Raw expectation of IY: 0.0026\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of YZ: -0.003\n",
      "Raw expectation of YI: 0.0068\n",
      "Raw expectation of IZ: 0.0008\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of ZX: 0.0044\n",
      "Raw expectation of ZI: 0.0104\n",
      "Raw expectation of IX: -0.0028\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of ZY: 0.0008\n",
      "Raw expectation of ZI: 0.0086\n",
      "Raw expectation of IY: 0.0032\n",
      "Raw expectation of II: 1.0\n",
      "Raw expectation of ZZ: 1.0\n",
      "Raw expectation of ZI: 0.0016\n",
      "Raw expectation of IZ: 0.0048\n",
      "Raw expectation of II: 1.0\n",
      "Expectations for x_x POVM: [1.0, -0.0122, -0.0052, 1.0]\n"
     ]
    }
   ],
   "source": [
    "observables = {}\n",
    "for povm_name in povm_names:\n",
    "    observables[povm_name] = obtain_expectations_for_qst(qc, p, generate_pauli_strings_from_povm_name(povm_name))\n",
    "\n",
    "print(f'Expectations for {povm_names[0]} POVM: {observables[povm_names[0]]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need empirical distributions that correspond to selected POVMs in order to perform QST.\n",
    "In ideal case of projective measurements, probability distribution of a POVM can be calculated from a set of Pauli observables and coefficient matrix $A_{\\mathrm{coefficient}}$.\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix} p_1\\\\ p_2 \\\\ \\vdots \\\\ p_{N}\\end{pmatrix}= A_{\\mathrm{coefficient}} \\begin{pmatrix} \\langle XX \\cdots X \\rangle \\\\ \\langle XX \\cdots I \\rangle \\\\  \\vdots \\\\ \\langle II \\cdots I \\rangle \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "Here we create a set of empirical distributions `empi_dists`, which are pairs of a repetition number and relative frequencies.\n",
    "We calculate `empi_dists` from measured expectations by using `calc_empi_dist_from_observables()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x_x: (10000, array([ 0.49565, -0.00175,  0.00175,  0.50435]))\n",
      "x_y: (10000, array([0.25625, 0.24645, 0.24385, 0.25345]))\n",
      "x_z: (10000, array([0.2536, 0.2438, 0.2453, 0.2573]))\n",
      "y_x: (10000, array([0.2476, 0.2444, 0.251 , 0.257 ]))\n",
      "y_y: (10000, array([ 0.00495,  0.50365,  0.49635, -0.00495]))\n",
      "y_z: (10000, array([0.25115, 0.25225, 0.24925, 0.24735]))\n",
      "z_x: (10000, array([0.253 , 0.2522, 0.2456, 0.2492]))\n",
      "z_y: (10000, array([0.25315, 0.25115, 0.24845, 0.24725]))\n",
      "z_z: (10000, array([ 0.5016, -0.0008,  0.0008,  0.4984]))\n"
     ]
    }
   ],
   "source": [
    "empi_dists = []\n",
    "for povm_name in povm_names:\n",
    "    empi_dist = calc_empi_dist_from_observables(observables[povm_name], num_shots, pauli_strings[povm_name], povms[povm_name])\n",
    "    empi_dists.append(empi_dist)\n",
    "    print(f\"{povm_name}: {empi_dist}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We generate a StandardQst object. From now on, the flow of QST will be the same as explained in [tutorial_for_standardqtomography](./tutorial_for_standardqtomography.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "tester_povms = []\n",
    "for povm_name in povm_names:\n",
    "    tester_povms.append(povms[povm_name])\n",
    "\n",
    "qst = StandardQst(tester_povms, on_para_eq_constraint=True, schedules=\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we choose a linear estimator, the estimate is calculated as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type:\n",
      "State\n",
      "\n",
      "Dim:\n",
      "4\n",
      "\n",
      "Vec:\n",
      "[ 5.00000000e-01 -1.80000000e-03  1.00000000e-03  5.66666667e-04\n",
      " -2.00000000e-03  5.00000000e-01  9.70000000e-03  1.09000000e-02\n",
      "  1.33333333e-03  4.60000000e-03 -5.00000000e-01 -1.50000000e-03\n",
      "  3.43333333e-03  2.20000000e-03  4.00000000e-04  5.00000000e-01]\n",
      "is estimate physical? :  False\n",
      "\n",
      "Eigenvalues are:  [1.0000779004210014, 0.007635233277463086, -0.0007271297383294317, -0.006986003960135253]\n"
     ]
    }
   ],
   "source": [
    "estimator = LinearEstimator()\n",
    "result = estimator.calc_estimate(qtomography=qst, empi_dists=empi_dists, is_computation_time_required=True)\n",
    "estimate = result.estimated_qoperation\n",
    "print(estimate)\n",
    "print(\"is estimate physical? : \", estimate.is_physical())\n",
    "print(\"\\nEigenvalues are: \", estimate.calc_eigenvalues())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An eigenvalue of the estimated density matrix is negative, which violates the requirement of positive-semidefiniteness on density matrix. This kind of violation can occur when we choose a linear estimator. In order to avoid the problem, we need to perform a constraint optimization at the data-processing."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "cb92e15a8d637e9dd9c5ecfc8cd130dbba6859f1dbeed3d128cdd82d23a19cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.8.7 64-bit ('venv': venv)",
   "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.8.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
