{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using Quara's standard tomography features from Forest (1-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": [
    "import numpy as np\n",
    "from numpy import pi\n",
    "\n",
    "# quara\n",
    "from quara.interface.forest.api import (\n",
    "    generate_preprocess_program,\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.objects.state_typical import generate_state_from_name\n",
    "from quara.protocol.qtomography.standard.standard_qst import StandardQst\n",
    "from quara.protocol.qtomography.standard.standard_qpt import StandardQpt\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, PHASE\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": [
    "## Quantum State Tomography (QST)"
   ]
  },
  {
   "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 1 qubits system\n",
    "qc = get_qc(\"1q-qvm\")\n",
    "qubits = [0]\n",
    "\n",
    "# define initialization of the quantum system\n",
    "# creating a state called \"a state\"\n",
    "num_shots = 10000\n",
    "p = Program()\n",
    "p += H(qubits[0])\n",
    "p += PHASE(pi/4, qubits[0])\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\n",
      "Type:\n",
      "Povm\n",
      "\n",
      "Dim:\n",
      "2\n",
      "\n",
      "Number of outcomes:\n",
      "2\n",
      "\n",
      "Vecs:\n",
      "[[ 0.70710678  0.70710678  0.          0.        ]\n",
      " [ 0.70710678 -0.70710678  0.          0.        ]]\n",
      "\n",
      "Pauli strings that will be used in this tomographical experiment:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'x': ['X', 'I'], 'y': ['Y', 'I'], 'z': ['Z', 'I']}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_sys = generate_composite_system(\"qubit\", 1)\n",
    "povms = {}\n",
    "pauli_strings = {}\n",
    "povm_names = [\"x\", \"y\", \"z\"]\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 X: 0.6976\n",
      "Raw expectation of I: 1.0\n",
      "Raw expectation of Y: 0.6948\n",
      "Raw expectation of I: 1.0\n",
      "Raw expectation of Z: 0.0024\n",
      "Raw expectation of I: 1.0\n",
      "Expectations for x POVM: [0.6976, 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: (10000, array([0.8488, 0.1512]))\n",
      "y: (10000, array([0.8474, 0.1526]))\n",
      "z: (10000, array([0.5012, 0.4988]))\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",
      "2\n",
      "\n",
      "Vec:\n",
      "[0.70710678 0.49327769 0.49129779 0.00169706]\n",
      "is estimate physical? :  True\n",
      "\n",
      "Eigenvalues are:  [0.9922901989680475, 0.007709801031952318]\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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Process Tomography"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, we consider process tomography on 1-qubit system.\n",
    "First, we define a `target_program` which is the quantum process we want to investigate. We also define a function `obtain_expectations_for_qpt` which is the function that calls pyquil backend to obtain observations. \n",
    "\n",
    "When performing process tomography, we need a pre process for state preparation which will be defined as `preprocess_program`. After executing `preprocess_program`, we run the `target_program` and then we measure the expectations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# example for 1 qubits system\n",
    "qc = get_qc(\"1q-qvm\")\n",
    "qubits = [0]\n",
    "num_shots = 10000\n",
    "\n",
    "# define a process that we want to estimate\n",
    "target_program = Program()\n",
    "target_program += PHASE(pi/4, qubits[0])\n",
    "\n",
    "def obtain_expectations_for_qpt(qc, num_shots, preprocess_program, 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_program = preprocess_program + program\n",
    "    tomo_experiment_program.wrap_in_numshots_loop(num_shots)\n",
    "    tomo_experiment = Experiment(settings, tomo_experiment_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": [
    "In quara, the name of a quantum system is expressed using strings such as `x0` or `z0_y0` in multiple qubit systems. There are 6 basic quantum state that are widely used in 1 qubit system which are the following:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "|\\psi_{z_0}\\rangle &= |0\\rangle \\\\\n",
    "|\\psi_{z_1}\\rangle &= |1\\rangle \\\\\n",
    "|\\psi_{x_0}\\rangle &= |0\\rangle + |1\\rangle \\\\\n",
    "|\\psi_{x_1}\\rangle &= |0\\rangle - |1\\rangle \\\\\n",
    "|\\psi_{y_0}\\rangle &= |0\\rangle + i|1\\rangle \\\\\n",
    "|\\psi_{y_1}\\rangle &= |0\\rangle - i|1\\rangle \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "We create a state object in Quara and then generate programs to realize those states using `generate_preprocess_program`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample of Quara tester object for state_name: x0\n",
      "Type:\n",
      "State\n",
      "\n",
      "Dim:\n",
      "2\n",
      "\n",
      "Vec:\n",
      "[0.70710678 0.70710678 0.         0.        ]\n",
      "\n",
      "Programs that will be used for state preparation in this tomographical experiment:\n",
      "x0: H 0\n",
      "\n",
      "y0: RX(-pi/2) 0\n",
      "\n",
      "z0: \n",
      "z1: X 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "c_sys = generate_composite_system(\"qubit\", 1)\n",
    "# state names of tester states\n",
    "state_names = [\"x0\", \"y0\", \"z0\", \"z1\"]\n",
    "states = {}\n",
    "preprocess_programs = {}\n",
    "for state_name in state_names:\n",
    "    states[state_name] = generate_state_from_name(c_sys, state_name)\n",
    "    # create pre-process programs for pyquil that correspond to quara's state\n",
    "    preprocess_programs[state_name] = generate_preprocess_program(qubits, state_name)\n",
    "print(f\"Sample of Quara tester object for state_name: {state_names[0]}\")\n",
    "print(states[state_names[0]])\n",
    "print(\"\\nPrograms that will be used for state preparation in this tomographical experiment:\")\n",
    "for state_name in state_names:\n",
    "    print(f\"{state_name}: {preprocess_programs[state_name]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alike quantum state tomography we generate Pauli strings and POVM for measurement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sample of Quara tester object for povm_name: x\n",
      "Type:\n",
      "Povm\n",
      "\n",
      "Dim:\n",
      "2\n",
      "\n",
      "Number of outcomes:\n",
      "2\n",
      "\n",
      "Vecs:\n",
      "[[ 0.70710678  0.70710678  0.          0.        ]\n",
      " [ 0.70710678 -0.70710678  0.          0.        ]]\n",
      "\n",
      "Pauli strings that will be used in this tomographical experiment:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'x': ['X', 'I'], 'y': ['Y', 'I'], 'z': ['Z', 'I']}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_sys = generate_composite_system(\"qubit\", 1)\n",
    "povms = {}\n",
    "pauli_strings = {}\n",
    "povm_names = [\"x\", \"y\", \"z\"]\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 obtain expectations for each POVM measurement and preprocess states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x0\n",
      "x\n",
      "Raw expectation of X: 0.7016\n",
      "Raw expectation of I: 1.0\n",
      "y\n",
      "Raw expectation of Y: 0.7198\n",
      "Raw expectation of I: 1.0\n",
      "z\n",
      "Raw expectation of Z: 0.0102\n",
      "Raw expectation of I: 1.0\n",
      "y0\n",
      "x\n",
      "Raw expectation of X: -0.7178\n",
      "Raw expectation of I: 1.0\n",
      "y\n",
      "Raw expectation of Y: 0.7034\n",
      "Raw expectation of I: 1.0\n",
      "z\n",
      "Raw expectation of Z: -0.0034\n",
      "Raw expectation of I: 1.0\n",
      "z0\n",
      "x\n",
      "Raw expectation of X: 0.0078\n",
      "Raw expectation of I: 1.0\n",
      "y\n",
      "Raw expectation of Y: -0.0018\n",
      "Raw expectation of I: 1.0\n",
      "z\n",
      "Raw expectation of Z: 1.0\n",
      "Raw expectation of I: 1.0\n",
      "z1\n",
      "x\n",
      "Raw expectation of X: -0.0012\n",
      "Raw expectation of I: 1.0\n",
      "y\n",
      "Raw expectation of Y: 0.0092\n",
      "Raw expectation of I: 1.0\n",
      "z\n",
      "Raw expectation of Z: -1.0\n",
      "Raw expectation of I: 1.0\n"
     ]
    }
   ],
   "source": [
    "observables = {}\n",
    "for state_name in state_names:\n",
    "    observables[state_name] = {}\n",
    "    print(state_name)\n",
    "    for povm_name in povm_names:\n",
    "        print(povm_name)\n",
    "        observable = obtain_expectations_for_qpt(qc, num_shots, preprocess_programs[state_name], target_program, pauli_strings[povm_name])\n",
    "        observables[state_name][povm_name] = observable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alike state tomography we calculate the empirical distribution using `calc_empi_dist_from_observation`.\n",
    "Make sure the orders of `empi_dists, tester_states, tester_povms` match to each other, otherwise the estimation cannot be caculated properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "empirical distributions for state preparation x0\n",
      "x: (10000, array([0.8508, 0.1492]))\n",
      "y: (10000, array([0.8599, 0.1401]))\n",
      "z: (10000, array([0.5051, 0.4949]))\n",
      "empirical distributions for state preparation y0\n",
      "x: (10000, array([0.1411, 0.8589]))\n",
      "y: (10000, array([0.8517, 0.1483]))\n",
      "z: (10000, array([0.4983, 0.5017]))\n",
      "empirical distributions for state preparation z0\n",
      "x: (10000, array([0.5039, 0.4961]))\n",
      "y: (10000, array([0.4991, 0.5009]))\n",
      "z: (10000, array([1., 0.]))\n",
      "empirical distributions for state preparation z1\n",
      "x: (10000, array([0.4994, 0.5006]))\n",
      "y: (10000, array([0.5046, 0.4954]))\n",
      "z: (10000, array([1.11022302e-16, 1.00000000e+00]))\n"
     ]
    }
   ],
   "source": [
    "empi_dists = []\n",
    "for state_name in state_names:\n",
    "    print(f\"empirical distributions for state preparation {state_name}\")\n",
    "    for povm_name in povm_names:\n",
    "        empi_dist = calc_empi_dist_from_observables(observables[state_name][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",
    "\n",
    "tester_states = []\n",
    "for state_name in state_names:\n",
    "    tester_states.append(states[state_name])\n",
    "\n",
    "tester_povms = []\n",
    "for povm_name in povm_names:\n",
    "    tester_povms.append(povms[povm_name])\n",
    "\n",
    "qpt = StandardQpt(states=tester_states, povms=tester_povms, on_para_eq_constraint=True, schedules=\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate is calculated as follows when linear estimator is chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type:\n",
      "Gate\n",
      "\n",
      "Dim:\n",
      "2\n",
      "\n",
      "HS:\n",
      "[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]\n",
      " [ 3.30000000e-03  6.98300000e-01 -7.21100000e-01  4.50000000e-03]\n",
      " [ 3.70000000e-03  7.16100000e-01  6.99700000e-01 -5.50000000e-03]\n",
      " [-1.11022302e-16  1.02000000e-02 -3.40000000e-03  1.00000000e+00]]\n",
      "is estimate physical? :  False\n",
      "\n",
      "Eigenvalues are: [-8.42148956e-03 -1.64624349e-03  7.55747516e-03  2.00251026e+00]\n"
     ]
    }
   ],
   "source": [
    "estimator = LinearEstimator()\n",
    "result = estimator.calc_estimate(qtomography=qpt, 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",
    "evals, evecs = np.linalg.eigh(estimate.to_choi_matrix())\n",
    "print(\"\\nEigenvalues are:\", evals)"
   ]
  }
 ],
 "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
}