{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An example of how to run GST on a 2-qubit system\n",
    "This example gives an overview of the typical steps used to perform an end-to-end (i.e. experimental-data-to-report) Gate Set Tomography analysis on a 2-qubit system.  The steps are very similar to the single-qubit case described in the tutorials, but we thought 2Q-GST is an important enough topic to deserve a separate example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import pygsti"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Construct the desired 2-qubit model\n",
    "Since the purpose of this example is to show how to *run* 2Q-GST, we'll just use a built-in \"standard\" 2-qubit model.  (Another example covers how to create a custom 2-qubit model.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pygsti.modelpacks import smq2Q_XYICNOT\n",
    "target_model = smq2Q_XYICNOT.target_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: create an experiment design\n",
    "An experiment design is a object containing all the information needed to perform and later interpret the data from a set of circuits.  In the case of GST, lists of fiducial and germ sub-circuits are the building blocks of the circuits performed in the experiment. Typically, these lists are either provided by pyGSTi because you're using a \"standard\" model (as we are here), or computed using the \"fiducial selection\" and \"germ selection\" algorithms which are a part of pyGSTi and covered in the tutorials.  As an additional input, we'll need a list of lengths indicating the maximum length circuits to use on each successive GST iteration.  Since 2Q-GST can take a while, only use short sequences (`max_max_lengths=1`) with fiducial-pair reduction (`fpr=True`) to demonstrate 2Q-GST more quickly (because we know you have important stuff to do)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_design = smq2Q_XYICNOT.get_gst_experiment_design(max_max_length=2, fpr=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3: Data generation\n",
    "Now that we have an experment design we can generate the list of experiments needed to run GST, just like in the 1-qubit case.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create an empty dataset file at example_files/My2QExample/data/dataset.txt, which stores the\n",
    "# list of experiments and zerod-out columns where data should be inserted.\n",
    "pygsti.io.write_empty_protocol_data(exp_design, \"example_files/My2QExample\", clobber_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate some \"fake\" (simulated) data based on a depolarized version of the target model.  In actual\n",
    "# situations, you'd fill in dataset.txt with real data.\n",
    "mdl_datagen = target_model.depolarize(op_noise=0.1, spam_noise=0.001)\n",
    "pygsti.io.fill_in_empty_dataset_with_fake_data(mdl_datagen, \"example_files/My2QExample/data/dataset.txt\",\n",
    "                                              nSamples=1000, seed=2020)\n",
    "\n",
    "# ---- NOTE: you can stop and restart the python session at this point; everything you need is saved to disk ---- \n",
    "\n",
    "#Load in the \"data object\" which packages together the dataset and experiment design\n",
    "data = pygsti.io.load_data_from_dir(\"example_files/My2QExample\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4: Run GST\n",
    "Just like for 1-qubit GST, we use the `StandardGST` protocol to compute the GST estimates.  Usually for two qubits this could take a long time (hours on a single cpu) based on the number of operation sequences used, and running on multiple processors is a good idea (see the MPI example).  Here, we set the tolerance to a high value ($10^{-3}$) so that it only takes around 30 minutes to run.\n",
    "\n",
    "Some notes about the options/arguments here that are particularly relevant to 2-qubit GST:\n",
    "  - `memlimit` gives an estimate of how much memory is available to use on your system (in bytes).  This is currently *not* a hard limit, and pyGSTi may require slightly more memory than this \"limit\".  So you'll need to be conservative in the value you place here: if your machine has 10GB of RAM, set this to 6 or 8 GB initially and increase it as you see how much memory is actually used using a separate OS performance monitor tool.  If you're running on multiple processors, this should be the memory available *per processor*.\n",
    "  - `verbosity` tells the routine how much detail to print to stdout.  If you don't mind waiting a while without getting any output, you can leave this at its default value (2).  If you can't standing wondering whether GST is still running or has locked up, set this to 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-- Std Practice:  Iter 1 of 1  (CPTP) --: \n",
      "  --- Iterative MLGST: Iter 1 of 2  907 operation sequences ---: \n",
      "    --- Minimum Chi^2 GST ---\n",
      "    Memory limit = 3.00GB\n",
      "    Cur, Persist, Gather = 0.15, 0.05, 0.29 GB\n",
      "      Evaltree generation (default) w/mem limit = 2.50GB\n",
      "      bulk_evaltree: created initial tree (907 strs) in 0s\n",
      "      bulk_evaltree: split tree (1 subtrees) in 0s\n",
      "       mem(1 subtrees, 1,1 param-grps, 1 proc-grps) in 0s = 3.33GB (3.33GB fc)\n",
      "      Created evaluation tree with 1 subtrees.  Will divide 1 procs into 1 (subtree-processing)\n",
      "       groups of ~1 procs each, to distribute over 1920 params (taken as 2 param groups of ~960 params).\n",
      "       Memory estimate = 1.66GB (cache=907, wrtLen1=960, wrtLen2=1920, subsPerProc=1).\n",
      "      --- Outer Iter 0: norm_f = 1.04145e+08, mu=1, |x|=1.09545e-07, |J|=4883.6\n",
      "      --- Outer Iter 1: norm_f = 192633, mu=305.73, |x|=1.3962, |J|=3116.42\n",
      "      --- Outer Iter 2: norm_f = 44469.9, mu=6522.24, |x|=0.936245, |J|=4420.69\n",
      "      --- Outer Iter 3: norm_f = 17253.1, mu=6093.18, |x|=0.830009, |J|=4265\n",
      "      --- Outer Iter 4: norm_f = 3216.4, mu=2031.06, |x|=0.77599, |J|=4735.84\n",
      "      --- Outer Iter 5: norm_f = 2654.34, mu=3597.57, |x|=0.798844, |J|=4783.64\n",
      "      --- Outer Iter 6: norm_f = 2317.62, mu=2646.96, |x|=0.818706, |J|=4850.27\n",
      "      --- Outer Iter 7: norm_f = 2146.59, mu=2220.23, |x|=0.841305, |J|=4900.44\n",
      "      --- Outer Iter 8: norm_f = 2036.34, mu=1908.06, |x|=0.858589, |J|=4950.73\n",
      "      --- Outer Iter 9: norm_f = 1965.23, mu=1657.54, |x|=0.870776, |J|=4998.76\n",
      "      --- Outer Iter 10: norm_f = 1918.97, mu=1459.25, |x|=0.880077, |J|=5039.71\n",
      "      --- Outer Iter 11: norm_f = 1887.87, mu=1314.15, |x|=0.887098, |J|=5073.96\n",
      "      --- Outer Iter 12: norm_f = 1865.34, mu=1179.47, |x|=0.892076, |J|=5103.07\n",
      "      --- Outer Iter 13: norm_f = 1848.35, mu=1041.57, |x|=0.895834, |J|=5127.92\n",
      "      --- Outer Iter 14: norm_f = 1834.91, mu=904.583, |x|=0.899, |J|=5149.42\n",
      "      --- Outer Iter 15: norm_f = 1823.7, mu=795.029, |x|=0.901841, |J|=5168.47\n",
      "      --- Outer Iter 16: norm_f = 1814.04, mu=709.82, |x|=0.904155, |J|=5185.93\n",
      "      --- Outer Iter 17: norm_f = 1806.25, mu=664.518, |x|=0.905889, |J|=5201.24\n",
      "      --- Outer Iter 18: norm_f = 1803.38, mu=708.108, |x|=0.907354, |J|=5210.24\n",
      "      --- Outer Iter 19: norm_f = 1797.11, mu=1368.31, |x|=0.9069, |J|=5221.67\n",
      "      --- Outer Iter 20: norm_f = 1793.14, mu=970.607, |x|=0.906877, |J|=5230.83\n",
      "      --- Outer Iter 21: norm_f = 1791.53, mu=970.733, |x|=0.907815, |J|=5234.12\n",
      "      --- Outer Iter 22: norm_f = 1789.32, mu=1674.4, |x|=0.907723, |J|=5240.32\n",
      "      Least squares message = Both actual and predicted relative reductions in the sum of squares are at most 0.001\n",
      "    Finding num_nongauge_params is too expensive: using total params.\n",
      "    Sum of Chi^2 = 1789.32 (2721 data params - 1920 model params = expected mean of 801; p-value = 0)\n",
      "    Completed in 1381.5s\n",
      "    2*Delta(log(L)) = 1796.15\n",
      "    Iteration 1 took 1381.7s\n",
      "    \n",
      "  --- Iterative MLGST: Iter 2 of 2  1082 operation sequences ---: \n",
      "    --- Minimum Chi^2 GST ---\n",
      "    Memory limit = 3.00GB\n",
      "    Cur, Persist, Gather = 0.28, 0.06, 0.29 GB\n",
      "      Evaltree generation (default) w/mem limit = 2.36GB\n",
      "      bulk_evaltree: created initial tree (1082 strs) in 1s\n",
      "      bulk_evaltree: split tree (1 subtrees) in 0s\n",
      "       mem(1 subtrees, 1,1 param-grps, 1 proc-grps) in 1s = 3.97GB (3.97GB fc)\n",
      "      Created evaluation tree with 1 subtrees.  Will divide 1 procs into 1 (subtree-processing)\n",
      "       groups of ~1 procs each, to distribute over 1920 params (taken as 2 param groups of ~960 params).\n",
      "       Memory estimate = 1.99GB (cache=1082, wrtLen1=960, wrtLen2=1920, subsPerProc=1).\n",
      "      --- Outer Iter 0: norm_f = 2553.12, mu=1, |x|=0.907723, |J|=5767.92\n",
      "      --- Outer Iter 1: norm_f = 2365.52, mu=2343.41, |x|=0.902086, |J|=5740.17\n",
      "      --- Outer Iter 2: norm_f = 2298.64, mu=1024.03, |x|=0.896497, |J|=5757.89\n",
      "      --- Outer Iter 3: norm_f = 2284.09, mu=1022.58, |x|=0.896905, |J|=5750.27\n",
      "      --- Outer Iter 4: norm_f = 2268.12, mu=556.574, |x|=0.895402, |J|=5762.38\n",
      "      --- Outer Iter 5: norm_f = 2265.56, mu=583.315, |x|=0.897171, |J|=5759.47\n",
      "      --- Outer Iter 6: norm_f = 2259.55, mu=534.095, |x|=0.896907, |J|=5767.02\n",
      "      --- Outer Iter 7: norm_f = 2257.27, mu=1060.51, |x|=0.896666, |J|=5770.45\n",
      "      --- Outer Iter 8: norm_f = 2256.01, mu=1060.4, |x|=0.896896, |J|=5772.09\n",
      "      Least squares message = Both actual and predicted relative reductions in the sum of squares are at most 0.001\n",
      "    Finding num_nongauge_params is too expensive: using total params.\n",
      "    Sum of Chi^2 = 2256.01 (3246 data params - 1920 model params = expected mean of 1326; p-value = 0)\n",
      "    Completed in 490.6s\n",
      "    2*Delta(log(L)) = 2263.46\n",
      "    Iteration 2 took 490.9s\n",
      "    \n",
      "    Switching to ML objective (last iteration)\n",
      "    --- MLGST ---\n",
      "    Memory: limit = 3.00GB(cur, persist, gthr = 0.36, 0.06, 0.29 GB)\n",
      "      --- Outer Iter 0: norm_f = 1131.73, mu=1, |x|=0.896896, |J|=4081.22\n",
      "      --- Outer Iter 1: norm_f = 1129.76, mu=754.219, |x|=0.89627, |J|=4085.7\n",
      "      Least squares message = Both actual and predicted relative reductions in the sum of squares are at most 0.001\n",
      "    Finding num_nongauge_params is too expensive: using total params.\n",
      "      Maximum log(L) = 1129.76 below upper bound of -2.41405e+06\n",
      "        2*Delta(log(L)) = 2259.53 (3246 data params - 1920 model params = expected mean of 1326; p-value = 0)\n",
      "      Completed in 101.7s\n",
      "    2*Delta(log(L)) = 2259.53\n",
      "    Final MLGST took 101.7s\n",
      "    \n",
      "  Iterative MLGST Total Time: 1974.3s\n",
      "      -- Performing 'stdgaugeopt' gauge optimization on CPTP estimate --\n",
      "Total time=0.570564 hours\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "start = time.time()\n",
    "protocol = pygsti.protocols.StandardGST(\"CPTP\", advancedOptions={'all': {'tolerance': 1e-3}}, verbosity=4)\n",
    "results = protocol.run(data, memlimit=3*(1024)**3)\n",
    "end = time.time()\n",
    "print(\"Total time=%f hours\" % ((end - start) / 3600.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5: Create report(s) using the returned `ModelEstimateResults` object\n",
    "The `ModelEstimateResults` object returned from `run` can be used to generate a \"general\" HTML report, just as in the 1-qubit case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running idle tomography\n",
      "Computing switchable properties\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/enielse/pyGSTi/pygsti/report/plothelpers.py:410: UserWarning:\n",
      "\n",
      "Number of degrees of freedom different for different boxes!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "report = pygsti.report.construct_standard_report(\n",
    "    results, title=\"Example 2Q-GST Report\", verbosity=2)\n",
    "report.write_html('example_files/easy_2q_report', verbosity=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now open [example_files/easy_2q_report/main.html](example_files/easy_2q_report/main.html) to see the results.  You've run 2-qubit GST!\n",
    "\n",
    "You can save the `ModelEstimateResults` object to the same directory as the data and experiment design:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.write()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}