{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Procedural Error Bars\n",
    "\n",
    "One other way we can use the `pygsti.report.reportables` module described in the [ModelAnalysisMetrics tutorial](ModelAnalysisMetrics.ipynb) is to procedurally generate error bars for any quantity you want.\n",
    "\n",
    "First, let's simulate a noisy GST experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygsti\n",
    "from pygsti.modelpacks import smq1Q_XY\n",
    "from pygsti.report import reportables as rptbl, modelfunction as modelfn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_model = smq1Q_XY.target_model()\n",
    "\n",
    "L=128\n",
    "edesign = smq1Q_XY.create_gst_experiment_design(L)\n",
    "\n",
    "noisy_model = target_model.randomize_with_unitary(.1)\n",
    "noisy_model = noisy_model.depolarize(.05)\n",
    "\n",
    "N=64\n",
    "dataset = pygsti.data.simulate_data(noisy_model,edesign,N)\n",
    "\n",
    "\n",
    "gst_proto = pygsti.protocols.StandardGST(modes=['full TP','CPTPLND','Target'],verbosity=2)\n",
    "data = pygsti.protocols.ProtocolData(edesign,dataset)\n",
    "results = gst_proto.run(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compute error bars on the CPTP estimate, and then get a 95% confidence interval \"view\" from the `ConfidenceRegionFactory`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crfact = results.estimates['CPTPLND'].add_confidence_region_factory('stdgaugeopt', 'final')\n",
    "crfact.compute_hessian(comm=None, mem_limit=3.0*(1024.0)**3) #optionally use multiple processors & set memlimit\n",
    "crfact.project_hessian('intrinsic error')\n",
    "\n",
    "crf_view = results.estimates['CPTPLND'].confidence_region_factories['stdgaugeopt','final'].view(95)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can construct `pygsti.report.ModelFunction` objects that take a function which computes some observable from a model and the extracted view from above to compute error bars on that quantity of interest.\n",
    "\n",
    "One common thing to check is error bars on the process matrices. The `ModelFunction` in this case only needs to return the operation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_model = results.estimates['CPTPLND'].models['stdgaugeopt'].copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_op(model, lbl):\n",
    "    return model[lbl]\n",
    "get_op_modelfn = modelfn.modelfn_factory(get_op)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rptbl.evaluate(get_op_modelfn(final_model, (\"Gxpi2\", 0)), crf_view)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rptbl.evaluate(get_op_modelfn(final_model, (\"Gypi2\", 0)), crf_view)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we can also create model functions that perform more complicated actions, such as computing other reportables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that when creating ModelFunctions in this way, the model where you want the quantity evaluated must be the first argument\n",
    "def ddist(model, ideal_model, lbl, basis):\n",
    "    return rptbl.half_diamond_norm(model[lbl], ideal_model[lbl], basis)\n",
    "ddist_modelfn = modelfn.modelfn_factory(ddist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rptbl.evaluate(ddist_modelfn(final_model, target_model, (\"Gxpi2\", 0), 'pp'), crf_view)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rptbl.evaluate(ddist_modelfn(final_model, target_model, (\"Gypi2\", 0), 'pp'), crf_view)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pygsti",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}