{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to generate error bars for 2Q-GST\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygsti\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 directory, run the 2Q-GST example)\n",
    "results = pygsti.io.read_results_from_dir(\"example_files/My2QExample\", \"StandardGST\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Layout creation w/mem limit = 3.00GB\n",
      "MatrixLayout: 1 processors divided into 1 x 1 x 1 (= 1) grid along circuit and parameter directions.\n",
      "   1 atoms, parameter block size limits (30, 30)\n",
      "*** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
      "    More atom-processors than hosts: each host gets ~1 atom-processors\n",
      "    Atom-processors already occupy a single node, dividing atom-processor into 1 param-processors.\n",
      "*** Divided 1-host atom-processor (~1 procs) into 1 param-processing groups ***\n",
      "    Param-processors already occupy a single node, dividing param-processor into 1 param2-processors.\n",
      "*** Divided 1-host param-processor (~1 procs) into 1 param2-processing groups ***\n",
      "   Esimated memory required = 1.9GB\n",
      "rank 0: 8.10623e-06s: beginning atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 27.3823s: block 1/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 56.8049s: block 2/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 86.7264s: block 3/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 116.186s: block 4/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 147.012s: block 5/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 176.884s: block 6/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 207.426s: block 7/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 237.48s: block 8/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 265.242s: block 9/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 293.074s: block 10/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 321.699s: block 11/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 349.566s: block 12/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 376.808s: block 13/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 404.092s: block 14/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 430.981s: block 15/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 458.278s: block 16/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 486.687s: block 17/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 513.101s: block 18/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 540.426s: block 19/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 568.835s: block 20/4096, atom 1/1, atom-size (#circuits) = 3884\n",
      "rank 0: 596.736s: block 21/4096, atom 1/1, atom-size (#circuits) = 3884\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-b70a459f9a3e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     11\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\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'CPTP'\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'stdgaugeopt'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'final'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\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[0mmem_limit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m3.0\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1024.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m#optionally use multiple processors & set memlimit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/protocols/confidenceregionfactory.py\u001b[0m in \u001b[0;36mcompute_hessian\u001b[0;34m(self, comm, mem_limit, approximate)\u001b[0m\n\u001b[1;32m    311\u001b[0m             \u001b[0mhessian_fn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_tools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogl_approximate_hessian\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mapproximate\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    312\u001b[0m                 \u001b[0;32melse\u001b[0m \u001b[0m_tools\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlogl_hessian\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 313\u001b[0;31m             hessian = hessian_fn(model, dataset, circuit_list,\n\u001b[0m\u001b[1;32m    314\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    315\u001b[0m                                  \u001b[0mcomm\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcomm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmem_limit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmem_limit\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[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/tools/likelihoodfns.py\u001b[0m in \u001b[0;36mlogl_hessian\u001b[0;34m(model, dataset, circuits, min_prob_clip, prob_clip_interval, radius, poisson_picture, op_label_aliases, mdc_store, comm, mem_limit, verbosity)\u001b[0m\n\u001b[1;32m    421\u001b[0m                          \u001b[0mregularization\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m'prob_clip_interval'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mprob_clip_interval\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    422\u001b[0m                          op_label_aliases, comm, mem_limit, ('hessian',), (), mdc_store, verbosity)\n\u001b[0;32m--> 423\u001b[0;31m     \u001b[0mlocal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# negative b/c objective is deltaLogL = max_logl - logL\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    424\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlayout\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mallgather_local_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'epp'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlocal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    425\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/objectivefns/objectivefns.py\u001b[0m in \u001b[0;36mhessian\u001b[0;34m(self, paramvec)\u001b[0m\n\u001b[1;32m   4877\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mex\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Hessian is not implemented for penalty terms yet!\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4878\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mparamvec\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfrom_vector\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mparamvec\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4879\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_gather_hessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_construct_hessian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcounts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtotal_counts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprob_clip_interval\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4880\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4881\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_hessian_from_block\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhprobs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdprobs12\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprobs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcounts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtotal_counts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfreqs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresource_alloc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/objectivefns/objectivefns.py\u001b[0m in \u001b[0;36m_construct_hessian\u001b[0;34m(self, counts, total_counts, prob_clip_interval)\u001b[0m\n\u001b[1;32m   1572\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[0mslicetup_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1573\u001b[0m                 \u001b[0mblk_rank\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparam2_resource_alloc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomm_rank\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1574\u001b[0;31m                 for (slice1, slice2, hprobs, dprobs12) in self.model.sim._iter_atom_hprobs_by_rectangle(\n\u001b[0m\u001b[1;32m   1575\u001b[0m                         atom, slicetup_list, True, param2_resource_alloc):\n\u001b[1;32m   1576\u001b[0m                     \u001b[0mlocal_slice1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_slct\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshift\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mslice1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mglobal_param_slice\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstart\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# indices into atom_hessian\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/forwardsims/distforwardsim.py\u001b[0m in \u001b[0;36m_iter_atom_hprobs_by_rectangle\u001b[0;34m(self, atom, wrt_slices_list, return_dprobs_12, resource_alloc)\u001b[0m\n\u001b[1;32m    321\u001b[0m             \u001b[0;31m# Note: no need to index w/ [atom.element_slice,...] (compare with _iter_hprobs_by_rectangles)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    322\u001b[0m             \u001b[0;31m# since these arrays are already sized to this particular atom (not to all the host's atoms)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 323\u001b[0;31m             self._bulk_fill_hprobs_dprobs_atom(hprobs, dprobs1, dprobs2, atom,\n\u001b[0m\u001b[1;32m    324\u001b[0m                                                wrtSlice1, wrtSlice2, resource_alloc)\n\u001b[1;32m    325\u001b[0m             \u001b[0;31m#Note: we give resource_alloc as our local `resource_alloc` above because all the arrays\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/forwardsims/distforwardsim.py\u001b[0m in \u001b[0;36m_bulk_fill_hprobs_dprobs_atom\u001b[0;34m(self, array_to_fill, deriv1_array_to_fill, deriv2_array_to_fill, atom, param_slice1, param_slice2, resource_alloc)\u001b[0m\n\u001b[1;32m    257\u001b[0m                                             param_slice2, resource_alloc)\n\u001b[1;32m    258\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 259\u001b[0;31m         self._bulk_fill_hprobs_atom(array_to_fill, host_param_slice1, host_param_slice2, atom,\n\u001b[0m\u001b[1;32m    260\u001b[0m                                     param_slice1, param_slice2, resource_alloc)\n\u001b[1;32m    261\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/pyGSTi/pygsti/forwardsims/matrixforwardsim.py\u001b[0m in \u001b[0;36m_bulk_fill_hprobs_atom\u001b[0;34m(self, array_to_fill, dest_param_slice1, dest_param_slice2, layout_atom, param_slice1, param_slice2, resource_alloc)\u001b[0m\n\u001b[1;32m   1435\u001b[0m             self._compute_dproduct_cache(layout_atom.tree, prodCache, scaleCache,\n\u001b[1;32m   1436\u001b[0m                                          resource_alloc, param_slice2)  # computed on rank=0 only\n\u001b[0;32m-> 1437\u001b[0;31m         hProdCache = self._compute_hproduct_cache(layout_atom.tree, prodCache, dProdCache1,\n\u001b[0m\u001b[1;32m   1438\u001b[0m                                                   \u001b[0mdProdCache2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mscaleCache\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresource_alloc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1439\u001b[0m                                                   param_slice1, param_slice2)  # computed on rank=0 only\n",
      "\u001b[0;32m~/pyGSTi/pygsti/forwardsims/matrixforwardsim.py\u001b[0m in \u001b[0;36m_compute_hproduct_cache\u001b[0;34m(self, layout_atom_tree, prod_cache, d_prod_cache1, d_prod_cache2, scale_cache, resource_alloc, wrt_slice1, wrt_slice2)\u001b[0m\n\u001b[1;32m    991\u001b[0m             \u001b[0mdLdR_sym\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdLdRa\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0m_np\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mswapaxes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdLdRb\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;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    992\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 993\u001b[0;31m             \u001b[0mhProdCache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miDest\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[0mhL\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mR\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mdLdR_sym\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0m_np\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtranspose\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[0mhR\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    994\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    995\u001b[0m             \u001b[0mscale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscale_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miDest\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mscale_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miLeft\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mscale_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0miRight\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mdot\u001b[0;34m(*args, **kwargs)\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",
    "# Choose forward simulator and how computation should be split up. Here are a couple examples:\n",
    "#results.estimates['CPTP'].models['stdgaugeopt'].sim = pygsti.forwardsims.MapForwardSimulator(num_atoms=100)\n",
    "results.estimates['CPTP'].models['stdgaugeopt'].sim = pygsti.forwardsims.MatrixForwardSimulator(param_blk_sizes=(30,30))\n",
    "\n",
    "# initialize a factory for the 'go0' gauge optimization within the 'default' estimate\n",
    "crfact = results.estimates['CPTP'].add_confidence_region_factory('stdgaugeopt', 'final')\n",
    "crfact.compute_hessian(comm=comm, mem_limit=3.0*(1024.0)**3) #optionally use multiple processors & set memlimit\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 disk\n",
    "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}