{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Intermediate Measurements Tutorial\n",
    "This tutorial will demonstrate how perform tomography on gate sets which, in addition to normal gates, contain *quantum instruments*.  Quantum instruments are maps that act on a qubit state (density matrix) and produce a qubit state along with a classical outcome.  That is, instruments are maps from $\\mathcal{B}(\\mathcal{H})$, the space of density matrices, to $\\mathcal{B}(\\mathcal{H}) \\otimes K(n)$, where $K(n)$ is a classical space of $n$ elements.\n",
    "\n",
    "In pyGSTi, instruments are represented as collections of gates, one for each classical \"outcome\" of the instrument.  This tutorial will demonstrate how to add instruments to `GateSet` objects, compute probabilities using such `GateSet`s, and ultimately perform tomography on them.  We'll start with a few familiar imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import pygsti\n",
    "from pygsti.construction import std1Q_XYI as std\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Instrument construction\n",
    "Next, we'll add an instrument to our \"standard\" gate set - a 1-qubit gate set containing $I$, $X(\\pi/2)$, and $Y(\\pi/2)$ gates.  The ideal instrument will be named `\"Iz\"` (all instrument names must begin with `\"I\"`), and consist of perfect projectors onto the 0 and 1 states.  Instead of labelling the associated outcomes \"0\" and \"1\", which might me most logical, we'll name them \"p0\" and \"p1\" so it's easier to distinguish them from the final POVM outcomes which *are* labelled \"0\" and \"1\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "#Make a copy so we don't modify the original\n",
    "gs_target = std.gs_target.copy()\n",
    "\n",
    "#Create and add the ideal instrument\n",
    "E0 = gs_target.effects['0']\n",
    "E1 = gs_target.effects['1']\n",
    " # Alternate indexing that uses POVM label explicitly\n",
    " # E0 = gs_target['Mdefault']['0']  # 'Mdefault' = POVM label, '0' = effect label\n",
    " # E1 = gs_target['Mdefault']['1']\n",
    "Gmz_plus = np.dot(E0,E0.T) #note effect vectors are stored as column vectors\n",
    "Gmz_minus = np.dot(E1,E1.T)\n",
    "gs_target['Iz'] = pygsti.obj.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n",
    "\n",
    "#For later use, record the identity POVM vector\n",
    "povm_ident = gs_target.effects['0'] + gs_target.effects['1'] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In order to generate some simulated data later on, we'll now create a noisy version of `gs_target` by depolarizing the state preparation, gates, and POVM, and also rotating the basis that is measured by the instrument and POVM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "gs_noisy = gs_target.depolarize(gate_noise=0.01, spam_noise=0.01)\n",
    "gs_noisy.effects.depolarize(0.01)  #because above call only depolarizes the state prep, not the POVMs\n",
    "\n",
    "# add a rotation error to the POVM\n",
    "Uerr = pygsti.rotation_gate_mx([0,0.02,0])\n",
    "gs_noisy.effects['0'] = np.dot(gs_noisy.effects['0'].T,Uerr).T\n",
    "gs_noisy.effects['1'] = povm_ident - gs_noisy.effects['0']\n",
    "\n",
    "#Could also do this:\n",
    "#E0 = np.dot(gs_noisy['Mdefault']['0'].T,Uerr).T\n",
    "#E1 = povm_ident - E0\n",
    "#gs_noisy['Mdefault'] = pygsti.obj.UnconstrainedPOVM({'0': E0, '1': E1})\n",
    "\n",
    "# Use the same rotated effect vectors to \"rotate\" the instrument Iz too\n",
    "E0 = gs_noisy.effects['0']\n",
    "E1 = gs_noisy.effects['1']\n",
    "Gmz_plus = np.dot(E0,E0.T)\n",
    "Gmz_minus = np.dot(E1,E1.T)\n",
    "gs_noisy['Iz'] = pygsti.obj.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n",
    "\n",
    "#print(gs_noisy) #print the gate set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Generating probabilities \n",
    "Instrument labels (e.g. `\"Iz\"`) may be included within `GateString` objects, and `GateSet` objects are able to compute probabilities for them just like normal (non-instrument) gate sequences.  The difference is that probabilities are labeled by tuples of instrument and POVM outcomes - referred to as **\"outcome tuples\"** - one for each instrument and one for the final POVM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('p0', '0'): 0.5000000000000003,\n",
       " ('p0', '1'): 0.0,\n",
       " ('p1', '0'): 0.0,\n",
       " ('p1', '1'): 0.4999999999999999}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dict(gs_target.probs( pygsti.obj.GateString(('Gx','Iz')) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{('p0', 'p0', '0'): 0.5000000000000006,\n",
       " ('p0', 'p0', '1'): 0.0,\n",
       " ('p0', 'p1', '0'): 0.0,\n",
       " ('p0', 'p1', '1'): 0.5000000000000001,\n",
       " ('p1', 'p0', '0'): 0.0,\n",
       " ('p1', 'p0', '1'): 0.0,\n",
       " ('p1', 'p1', '0'): 0.0,\n",
       " ('p1', 'p1', '1'): 0.0}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dict(gs_target.probs( pygsti.obj.GateString(('Iz','Gx','Iz')) ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In fact, pyGSTi *always* labels probabilties using outcome tuples, it's just that in the non-instrument case they're always 1-tuples and by `OutcomeLabelDict` magic can be treated as if they were just strings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probs =  {('0',): 0.5000000000000002, ('1',): 0.4999999999999998}\n",
      "probs['0'] =  0.5000000000000002\n",
      "probs[('0',)] =  0.5000000000000002\n"
     ]
    }
   ],
   "source": [
    "probs = gs_target.probs( pygsti.obj.GateString(('Gx',)) )\n",
    "print(\"probs = \",dict(probs))\n",
    "print(\"probs['0'] = \", probs['0']) #This works...\n",
    "print(\"probs[('0',)] = \", probs[('0',)]) # and so does this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Performing tomography\n",
    "\n",
    "### Simulated data generation\n",
    "Now let's perform tomography on a gate set that includes instruments.  First, we'll generate some data using `gs_noisy` in exactly the same way as we would for any other gate set:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "germs = std.germs\n",
    "fiducials = std.fiducials\n",
    "max_lengths = [1] # keep it simple & fast\n",
    "\n",
    "lsgst_list = pygsti.construction.make_lsgst_experiment_list(\n",
    "    gs_noisy,fiducials,fiducials,germs,max_lengths)\n",
    "\n",
    "#print(\"Gate sequences:\")\n",
    "#print(lsgst_list) #note that this contains LGST strings with \"Iz\"\n",
    "\n",
    "#Create the DataSet\n",
    "ds = pygsti.construction.generate_fake_data(gs_noisy,lsgst_list,1000,'multinomial',seed=2018)\n",
    "\n",
    "#Write it to a text file to demonstrate the format:\n",
    "pygsti.io.write_dataset(\"tutorial_files/intermediate_meas_dataset.txt\",ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Notice the format of [tutorial_files/intermediate_meas_dataset.txt](tutorial_files/intermediate_meas_dataset.txt), which includes a column for each distinct outcome tuple.  Since not all experiments contain data for all outcome tuples, the `\"--\"` is used as a placeholder.  Now that the data is generated, we run LGST or LSGST just like we would for any other gateset:\n",
    "\n",
    "### LGST"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gateset Difference:\n",
      " Preps:\n",
      "  rho0 = 0.0293321\n",
      " POVMs:\n",
      "  Mdefault:     0 = 0.0199685\n",
      "    1 = 0.0184818\n",
      " Gates:\n",
      "  Gi = 0.0614996\n",
      "  Gx = 0.0397322\n",
      "  Gy = 0.0253245\n",
      " Instruments:\n",
      "  Iz:     p0 = 0.0393195\n",
      "    p1 = 0.0322578\n",
      "\n",
      "Frobdiff after GOpt =  0.010524022143914753\n"
     ]
    }
   ],
   "source": [
    "#Run LGST\n",
    "gs_lgst = pygsti.do_lgst(ds, fiducials,fiducials, gs_target)\n",
    "#print(gs_lgst)\n",
    "\n",
    "#Gauge optimize the result to the true data-generating gate set (gs_noisy),\n",
    "# and compare.  Mismatch is due to finite sample noise.\n",
    "gs_lgst_opt = pygsti.gaugeopt_to_target(gs_lgst,gs_noisy)\n",
    "print(gs_noisy.strdiff(gs_lgst_opt))\n",
    "print(\"Frobdiff after GOpt = \",gs_noisy.frobeniusdist(gs_lgst_opt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Long-sequence GST\n",
    "Instruments just add parameters to a `GateSet` like gates, state preparations, and POVMs do.  The total number of parameters in our gate set is \n",
    "\n",
    "$4$ (prep) + $2\\times 4$ (2 effects) + $5\\times 16$ (3 gates and 2 instrument members) $ = 92$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "92"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gs_target.num_params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Gate Sequence Creation ---\n",
      "   92 sequences created\n",
      "   Dataset has 128 entries: 92 utilized, 0 requested sequences were missing\n",
      "--- LGST ---\n",
      "  Singular values of I_tilde (truncating to first 4 of 6) = \n",
      "  4.243132186580752\n",
      "  1.3675281047825063\n",
      "  1.3451801562592671\n",
      "  1.320952449823206\n",
      "  0.052840411822390194\n",
      "  0.020469323956352266\n",
      "  \n",
      "  Singular values of target I_tilde (truncating to first 4 of 6) = \n",
      "  4.242640687119285\n",
      "  1.4142135623730954\n",
      "  1.4142135623730947\n",
      "  1.4142135623730945\n",
      "  3.1723744950054595e-16\n",
      "  1.0852733691121267e-16\n",
      "  \n",
      "--- Iterative MLGST: Iter 1 of 1  92 gate strings ---: \n",
      "  --- Minimum Chi^2 GST ---\n",
      "  Sum of Chi^2 = 56.2886 (92 data params - 76 model params = expected mean of 16; p-value = 2.18045e-06)\n",
      "  Completed in 0.2s\n",
      "  2*Delta(log(L)) = 56.0781\n",
      "  Iteration 1 took 0.2s\n",
      "  \n",
      "  Switching to ML objective (last iteration)\n",
      "  --- MLGST ---\n",
      "    Maximum log(L) = 27.9886 below upper bound of -138546\n",
      "      2*Delta(log(L)) = 55.9771 (92 data params - 76 model params = expected mean of 16; p-value = 2.45485e-06)\n",
      "    Completed in 0.1s\n",
      "  2*Delta(log(L)) = 55.9771\n",
      "  Final MLGST took 0.2s\n",
      "  \n",
      "Iterative MLGST Total Time: 0.4s\n",
      "  -- Adding Gauge Optimized (go0) --\n",
      "--- Re-optimizing logl after robust data scaling ---\n",
      "  --- MLGST ---\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/enielse/research/pyGSTi/packages/pygsti/objects/estimate.py:513: UserWarning:\n",
      "\n",
      "Max-model params (92) <= gate set params (92)!  Using k == 1.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    Maximum log(L) = 27.9886 below upper bound of -138546\n",
      "      2*Delta(log(L)) = 55.9771 (92 data params - 76 model params = expected mean of 16; p-value = 2.45485e-06)\n",
      "    Completed in 0.1s\n",
      "  -- Adding Gauge Optimized (go0) --\n"
     ]
    }
   ],
   "source": [
    "#Run long sequence GST\n",
    "results = pygsti.do_long_sequence_gst(ds,gs_target,fiducials,fiducials,germs,max_lengths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Frobdiff after GOpt =  0.008310755067937731\n"
     ]
    }
   ],
   "source": [
    "#Compare estimated gate set (after gauge opt) to data-generating one\n",
    "gs_est = results.estimates['default'].gatesets['go0']\n",
    "gs_est_opt = pygsti.gaugeopt_to_target(gs_est,gs_noisy)\n",
    "print(\"Frobdiff after GOpt = \", gs_noisy.frobeniusdist(gs_est_opt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The same analysis can be done for a trace-preserving gate set, whose instruments are constrained to *add* to a perfectly trace-preserving map.  The number of parameters in the gate set are now:  \n",
    "\n",
    "$3$ (prep) + $1\\times 4$ (effect and complement) + $3\\times 12$ (3 gates) + $(2\\times 16 - 3)$ (TP instrument) $ = 71$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POVM type =  <class 'pygsti.objects.povm.TPPOVM'>  Np= 4\n",
      "Instrument type =  <class 'pygsti.objects.instrument.TPInstrument'>  Np= 28\n",
      "Number of gateset parameters =  71\n"
     ]
    }
   ],
   "source": [
    "gs_targetTP = gs_target.copy()\n",
    "gs_targetTP.set_all_parameterizations(\"TP\")\n",
    "print(\"POVM type = \",type(gs_targetTP[\"Mdefault\"]),\" Np=\",gs_targetTP[\"Mdefault\"].num_params())\n",
    "print(\"Instrument type = \",type(gs_targetTP[\"Iz\"]),\" Np=\",gs_targetTP[\"Iz\"].num_params())\n",
    "print(\"Number of gateset parameters = \", gs_targetTP.num_params())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Gate Sequence Creation ---\n",
      "   92 sequences created\n",
      "   Dataset has 128 entries: 92 utilized, 0 requested sequences were missing\n",
      "--- LGST ---\n",
      "  Singular values of I_tilde (truncating to first 4 of 6) = \n",
      "  4.243132186580752\n",
      "  1.3675281047825063\n",
      "  1.3451801562592671\n",
      "  1.320952449823206\n",
      "  0.052840411822390194\n",
      "  0.020469323956352266\n",
      "  \n",
      "  Singular values of target I_tilde (truncating to first 4 of 6) = \n",
      "  4.242640687119285\n",
      "  1.4142135623730954\n",
      "  1.4142135623730947\n",
      "  1.4142135623730945\n",
      "  3.1723744950054595e-16\n",
      "  1.0852733691121267e-16\n",
      "  \n",
      "--- Iterative MLGST: Iter 1 of 1  92 gate strings ---: \n",
      "  --- Minimum Chi^2 GST ---\n",
      "  Sum of Chi^2 = 59.2126 (92 data params - 63 model params = expected mean of 29; p-value = 0.000773137)\n",
      "  Completed in 0.1s\n",
      "  2*Delta(log(L)) = 59.21\n",
      "  Iteration 1 took 0.2s\n",
      "  \n",
      "  Switching to ML objective (last iteration)\n",
      "  --- MLGST ---\n",
      "    Maximum log(L) = 29.555 below upper bound of -138546\n",
      "      2*Delta(log(L)) = 59.11 (92 data params - 63 model params = expected mean of 29; p-value = 0.00079598)\n",
      "    Completed in 0.1s\n",
      "  2*Delta(log(L)) = 59.11\n",
      "  Final MLGST took 0.2s\n",
      "  \n",
      "Iterative MLGST Total Time: 0.4s\n",
      "  -- Adding Gauge Optimized (go0) --\n",
      "--- Re-optimizing logl after robust data scaling ---\n",
      "  --- MLGST ---\n",
      "    Maximum log(L) = 29.555 below upper bound of -138546\n",
      "      2*Delta(log(L)) = 59.11 (92 data params - 63 model params = expected mean of 29; p-value = 0.00079598)\n",
      "    Completed in 0.1s\n",
      "  -- Adding Gauge Optimized (go0) --\n"
     ]
    }
   ],
   "source": [
    "resultsTP = pygsti.do_long_sequence_gst(ds,gs_targetTP,fiducials,fiducials,germs,max_lengths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Frobdiff after GOpt =  0.008137806725675554\n"
     ]
    }
   ],
   "source": [
    "#Again compare estimated gate set (after gauge opt) to data-generating one\n",
    "gs_est = resultsTP.estimates['default'].gatesets['go0']\n",
    "gs_est_opt = pygsti.gaugeopt_to_target(gs_est,gs_noisy)\n",
    "print(\"Frobdiff after GOpt = \", gs_noisy.frobeniusdist(gs_est_opt))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "**Thats it!**  You've done tomography on a gate set with intermediate measurments (instruments)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}