{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An example showing how to generate bootstrapped error bars."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import time\n",
    "import json\n",
    "import numpy as np\n",
    "import pygsti\n",
    "from pygsti.modelpacks import smq1Q_XY\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get a GST estimate (similar to Tutorial 0)\n",
    "\n",
    "# 1) get the target Model\n",
    "target_model = smq1Q_XY.target_model('full TP')\n",
    "\n",
    "# 2) get the building blocks needed to specify which operation sequences are needed\n",
    "prep_fiducials = smq1Q_XY.prep_fiducials()[0:4]\n",
    "meas_fiducials = smq1Q_XY.meas_fiducials()[0:3]\n",
    "germs = smq1Q_XY.germs()\n",
    "maxLengths = [1,2,4,8]\n",
    "\n",
    "# 3) generate \"fake\" data from a depolarized version of target_model\n",
    "mdl_datagen = target_model.depolarize(op_noise=0.1, spam_noise=0.001)\n",
    "listOfExperiments = pygsti.circuits.create_lsgst_circuits(\n",
    "    target_model, prep_fiducials, meas_fiducials, germs, maxLengths)\n",
    "ds = pygsti.data.simulate_data(mdl_datagen, listOfExperiments, num_samples=1000,\n",
    "                                            sample_error=\"binomial\", seed=1234)\n",
    "\n",
    "\n",
    "results = pygsti.run_stdpractice_gst(ds, target_model, prep_fiducials, meas_fiducials,\n",
    "                                    germs, maxLengths, modes=\"full TP\")\n",
    "estimated_model = results.estimates['full TP'].models['stdgaugeopt']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parametric Bootstrapping\n",
    "Here we do parametric bootstrapping, as indicated by the 'parametric' argument below.\n",
    "The output is eventually stored in the \"mean\" and \"std\" Models, which hold the mean and standard deviation values of the set of bootstrapped models (after gauge optimization).  It is this latter \"standard deviation Model\"\n",
    "which holds the collection of error bars.  Note: due to print setting issues, the outputs that are printed here will not necessarily reflect the true accuracy of the estimates made.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#The number of simulated datasets & models made for bootstrapping purposes.  \n",
    "# For good statistics, should probably be greater than 10.\n",
    "numGatesets=10\n",
    "\n",
    "param_boot_models = pygsti.drivers.create_bootstrap_models(\n",
    "                        numGatesets, ds, 'parametric', prep_fiducials, meas_fiducials, germs, maxLengths,\n",
    "                        input_model=estimated_model, target_model=target_model,\n",
    "                        start_seed=0, return_data=False,\n",
    "                        verbosity=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gauge_opt_pboot_models = pygsti.drivers.gauge_optimize_models(param_boot_models, estimated_model,\n",
    "                                                                 plot=False) #plotting support removed w/matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gauge_opt_pboot_models[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pboot_mean = pygsti.drivers.bootstrap._to_mean_model(gauge_opt_pboot_models, estimated_model)\n",
    "pboot_std  = pygsti.drivers.bootstrap._to_std_model(gauge_opt_pboot_models, estimated_model)\n",
    "\n",
    "#Summary of the error bars\n",
    "print(\"Parametric bootstrapped error bars, with\", numGatesets, \"resamples\\n\")\n",
    "print(\"Error in rho vec:\") \n",
    "print(pboot_std['rho0'].to_vector(), end='\\n\\n')\n",
    "print(\"Error in effect vecs:\")\n",
    "print(pboot_std['Mdefault'].to_vector(), end='\\n\\n')\n",
    "print(\"Error in Gxpi2:\")\n",
    "print(pboot_std['Gxpi2',0].to_vector(), end='\\n\\n')\n",
    "print(\"Error in Gypi2:\")\n",
    "print(pboot_std['Gypi2',0].to_vector())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-parametric Bootstrapping\n",
    "Here we do non-parametric bootstrapping, as indicated by the 'nonparametric' argument below.\n",
    "The output is again eventually stored in the \"mean\" and \"std\" Models, which hold the mean and standard deviation values of the set of bootstrapped models (after gauge optimization).  It is this latter \"standard deviation Model\"\n",
    "which holds the collection of error bars.  Note: due to print setting issues, the outputs that are printed here will not necessarily reflect the true accuracy of the estimates made.\n",
    "\n",
    "(Technical note: ddof = 1 is by default used when computing the standard deviation -- see numpy.std -- meaning that we are computing a standard deviation of the sample, not of the population.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#The number of simulated datasets & models made for bootstrapping purposes.  \n",
    "# For good statistics, should probably be greater than 10.\n",
    "numModels=10\n",
    "\n",
    "nonparam_boot_models = pygsti.drivers.create_bootstrap_models(\n",
    "                          numModels, ds, 'nonparametric', prep_fiducials, meas_fiducials, germs, maxLengths,\n",
    "                          target_model=target_model, start_seed=0, return_data=False, verbosity=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gauge_opt_npboot_models = pygsti.drivers.gauge_optimize_models(nonparam_boot_models, estimated_model,\n",
    "                                                                 plot=False) #plotting removed w/matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npboot_mean = pygsti.drivers.bootstrap._to_mean_model(gauge_opt_npboot_models, estimated_model)\n",
    "npboot_std  = pygsti.drivers.bootstrap._to_std_model(gauge_opt_npboot_models, estimated_model)\n",
    "\n",
    "#Summary of the error bars\n",
    "print(\"Non-parametric bootstrapped error bars, with\", numGatesets, \"resamples\\n\")\n",
    "print(\"Error in rho vec:\")\n",
    "print(npboot_std['rho0'].to_vector(), end='\\n\\n')\n",
    "print(\"Error in effect vecs:\")\n",
    "print(npboot_std['Mdefault'].to_vector(), end='\\n\\n')\n",
    "print(\"Error in Gxpi2:\")\n",
    "print(npboot_std['Gxpi2',0].to_vector(), end='\\n\\n')\n",
    "print(\"Error in Gypi2:\")\n",
    "print(npboot_std['Gypi2',0].to_vector())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.loglog(npboot_std.to_vector(),pboot_std.to_vector(),'.')\n",
    "plt.loglog(np.logspace(-4,-2,10),np.logspace(-4,-2,10),'--')\n",
    "plt.xlabel('Non-parametric')\n",
    "plt.ylabel('Parametric')\n",
    "plt.xlim((1e-4,1e-2)); plt.ylim((1e-4,1e-2))\n",
    "plt.title('Scatter plot comparing param vs. non-param bootstrapping error bars.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "api_updates",
   "language": "python",
   "name": "api_updates"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}