{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to generate error bars for 2Q-GST\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Memory limit was =  3221225472\n",
      "Memory limit is now =  2684354560.0\n"
     ]
    }
   ],
   "source": [
    "import pygsti\n",
    "import pickle\n",
    "import time\n",
    "\n",
    "#If we were using MPI\n",
    "# from mpi4py import MPI\n",
    "# comm = MPI.COMM_WORLD\n",
    "comm = None\n",
    "\n",
    "#Load the 2-qubit results (if you don't have this file, run the 2Q-GST example)\n",
    "with open(\"example_files/easy_2q_results.pkl\",\"rb\") as f:\n",
    "    results = pickle.load(f)\n",
    "\n",
    "#Set a memory limit\n",
    "print(\"Memory limit was = \", results.estimates['default'].parameters.get('memLimit',\"none given\"))\n",
    "results.estimates['default'].parameters['memLimit'] = 2.5*(1024.0)**3 # 2.5GB\n",
    "print(\"Memory limit is now = \", results.estimates['default'].parameters['memLimit'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaltree generation (deriv) w/mem limit = 2.48GB\n",
      " mem(1 subtrees, 1,1 param-grps, 1 proc-grps) in 0s = 6773.17GB (6773.17GB 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 (1616,1616) params (taken as 1616,4 param groups of ~1,404 params).\n",
      " Memory estimate = 2.08GB (cache=1317, wrtLen1=1, wrtLen2=404, subsPerProc=1).\n",
      "rank 0: 28.7456s: block 0/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 57.0849s: block 1/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 85.307s: block 2/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 114.001s: block 3/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 142.315s: block 4/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 170.56s: block 5/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 198.853s: block 6/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 227.573s: block 7/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 255.879s: block 8/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 284.05s: block 9/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 312.368s: block 10/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 340.616s: block 11/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 368.909s: block 12/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 397.244s: block 13/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 425.513s: block 14/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 453.775s: block 15/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 482.084s: block 16/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 510.291s: block 17/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 538.496s: block 18/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 566.854s: block 19/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 595.287s: block 20/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 623.537s: block 21/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 651.747s: block 22/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 679.983s: block 23/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 708.176s: block 24/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 736.482s: block 25/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 764.708s: block 26/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 792.868s: block 27/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 821.206s: block 28/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 849.412s: block 29/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 877.614s: block 30/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 905.894s: block 31/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 934.147s: block 32/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 962.373s: block 33/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 990.554s: block 34/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1018.69s: block 35/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1046.94s: block 36/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1075.31s: block 37/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1103.54s: block 38/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1131.84s: block 39/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1160.28s: block 40/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1188.52s: block 41/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1216.72s: block 42/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1244.92s: block 43/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1273.29s: block 44/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1301.41s: block 45/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1329.71s: block 46/4043, sub-tree 0/1, sub-tree-len = 1317\n",
      "rank 0: 1358.43s: block 47/4043, sub-tree 0/1, sub-tree-len = 1317\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-2-7208f561b6ac>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;31m# initialize a factory for the 'go0' gauge optimization within the 'default' estimate\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0mcrfact\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mresults\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mestimates\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'default'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_confidence_region_factory\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'go0'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'final'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0mcrfact\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompute_hessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcomm\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcomm\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m#optionally use multiple processors\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m \u001b[0mcrfact\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mproject_hessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'intrinsic error'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/research/pyGSTi/packages/pygsti/objects/confidenceregionfactory.py\u001b[0m in \u001b[0;36mcompute_hessian\u001b[0;34m(self, comm, memLimit, approximate)\u001b[0m\n\u001b[1;32m    239\u001b[0m                                  \u001b[0mminProbClip\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprobClipInterval\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mradius\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    240\u001b[0m                                  \u001b[0mcomm\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcomm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmemLimit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmemLimit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbosity\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mvb\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 241\u001b[0;31m                                  gateLabelAliases=aliases)\n\u001b[0m\u001b[1;32m    242\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    243\u001b[0m             nonMarkRadiusSq = max( 2*(_tools.logl_max(gateset, dataset)\n",
      "\u001b[0;32m~/research/pyGSTi/packages/pygsti/tools/likelihoodfns.py\u001b[0m in \u001b[0;36mlogl_hessian\u001b[0;34m(gateset, dataset, gatestring_list, minProbClip, probClipInterval, radius, poissonPicture, check, comm, memLimit, gateLabelAliases, smartc, verbosity)\u001b[0m\n\u001b[1;32m    733\u001b[0m         \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mkmax\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmySliceTupList\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    734\u001b[0m         for (slice1,slice2,hprobs,dprobs12) in gateset.bulk_hprobs_by_block(\n\u001b[0;32m--> 735\u001b[0;31m             evalSubTree, mySliceTupList, True, blkComm):\n\u001b[0m\u001b[1;32m    736\u001b[0m             \u001b[0mrank\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcomm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mGet_rank\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mcomm\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    737\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/research/pyGSTi/packages/pygsti/objects/gatematrixcalc.py\u001b[0m in \u001b[0;36mbulk_hprobs_by_block\u001b[0;34m(self, evalTree, wrtSlicesList, bReturnDProbs12, comm)\u001b[0m\n\u001b[1;32m   2683\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2684\u001b[0m                 dProdCache2 =self._compute_dproduct_cache(\n\u001b[0;32m-> 2685\u001b[0;31m                     evalTree, prodCache, scaleCache, comm, wrtSlice2)\n\u001b[0m\u001b[1;32m   2686\u001b[0m                 \u001b[0mdGs2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mevalTree\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfinal_view\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdProdCache2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2687\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/research/pyGSTi/packages/pygsti/objects/gatematrixcalc.py\u001b[0m in \u001b[0;36m_compute_dproduct_cache\u001b[0;34m(self, evalTree, prodCache, scaleCache, comm, wrtSlice, profiler)\u001b[0m\n\u001b[1;32m    962\u001b[0m             \u001b[0mL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mR\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprodCache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miLeft\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprodCache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miRight\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    963\u001b[0m             \u001b[0mdL\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdR\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdProdCache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miLeft\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdProdCache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miRight\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 964\u001b[0;31m             \u001b[0mdProdCache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_np\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mR\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    965\u001b[0m                 \u001b[0m_np\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mswapaxes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_np\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdR\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m#dot(dS, T) + dot(S, dT)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    966\u001b[0m             \u001b[0mprofiler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_time\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"compute_dproduct_cache: dots\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# error bars in reports require the presence of a fully-initialized\n",
    "# \"confidence region factory\" within the relevant Estimate object.\n",
    "# In most cases \"fully-initialized\" means that a Hessian has been \n",
    "# computed and projected onto the non-gauge space.\n",
    "start = time.time()\n",
    "\n",
    "# initialize a factory for the 'go0' gauge optimization within the 'default' estimate\n",
    "crfact = results.estimates['default'].add_confidence_region_factory('go0', 'final')\n",
    "crfact.compute_hessian(comm=comm) #optionally use multiple processors\n",
    "crfact.project_hessian('intrinsic error')\n",
    "\n",
    "end = time.time()\n",
    "print(\"Total time=%f hours\" % ((end - start) / 3600.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note above cell was executed for demonstration purposes, and was **keyboard-interrupted intentionally** since it would have taken forever on a single processor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#write results back to file\n",
    "with open(\"example_files/easy_2q_results_withCI.pkl\",\"wb\") as f:\n",
    "    pickle.dump(results, f)"
   ]
  },
  {
   "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": 2
}