Repository: sandialabs/pyGSTi (language: python)

Issue:
Evotype dependent serialization bugs
Some failing unit tests on the GST checkpointing branch revealed a serialization bug (really a bug that becomes apparent when trying to do serialization) that seems to be related to the use of the 'statevec' evotype. In particular, it looks like when using the 'statevec' evotype the _basis attribute is not always set properly remains initialized to None. We then encounter an error in serialization which expects there to be a basis object present and fails when trying to call the expected basis object's to_nice_serialization method.

Minimal working example for reproducing the bug:

```
from pygsti.modelpacks import smq1Q_XYI as std
test_model = std.target_model("static unitary", evotype='statevec')
test_model.to_nice_serialization()
```

Which returns the traceback:

```
 ---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 1
----> 1 test_model.to_nice_serialization()

File ~\Documents\pyGSTi_master_branch\pygsti\baseobjs\nicelyserializable.py:164, in NicelySerializable.to_nice_serialization(self)
    154 def to_nice_serialization(self):
    155     """
    156     Serialize this object in a way that adheres to "niceness" rules of common text file formats.
    157 
   (...)
    162         Python type, e.g. a string or list.
    163     """
--> 164     state = self._to_nice_serialization()
    165     if self._dbcoordinates is not None and isinstance(state, dict):
    166         # HACK: the code below adds the special dbcoordinates attribute to `state` when we're serializing to
    167         # something other than a MongoDB (which holds the dbcoordinates intrinsically in the document's
   (...)
    170         # a flag or something else that determines this, but for now we add this fragile HACK that checks the
    171         # caller names for 'add_mongodb_write_ops' (only in the call stack when we're serializing to MongoDB).
    172         import inspect

File ~\Documents\pyGSTi_master_branch\pygsti\models\explicitmodel.py:1603, in ExplicitOpModel._to_nice_serialization(self)
   1584 state.update({'basis': self.basis.to_nice_serialization(),
   1585               'default_gate_type': self.operations.default_param,
   1586               'default_prep_type': self.preps.default_param,
   (...)
   1599                                        if (self._param_interposer is not None) else None)
   1600               })
   1602 mmgraph = self.create_modelmember_graph()
-> 1603 state['modelmembers'] = mmgraph.create_serialization_dict()
   1604 return state

File ~\Documents\pyGSTi_master_branch\pygsti\modelmembers\modelmembergraph.py:181, in ModelMemberGraph.create_serialization_dict(self)
    179 sdict = OrderedDict()
    180 for mm_id, mm_node in self.mm_memo.items():
--> 181     sdict[str(mm_node.serialize_id)] = mm_node.mm.to_memoized_dict(self.mm_memo)
    183 # Tag extra metadata for root nodes
    184 for mm_type, root_nodes in self.mm_nodes.items():

File ~\Documents\pyGSTi_master_branch\pygsti\modelmembers\states\densestate.py:352, in DensePureState.to_memoized_dict(self, mmg_memo)
    349 mm_dict = super().to_memoized_dict(mmg_memo)
    351 mm_dict['dense_state_vector'] = self._encodemx(self.to_dense('Hilbert'))
--> 352 mm_dict['basis'] = self._basis.to_nice_serialization()
    354 return mm_dict

AttributeError: 'NoneType' object has no attribute 'to_nice_serialization'

```

A few additional notes. This looks to be specifically related to the use of the 'statevec' evotype. Excluding the evotype kwarg (I'm guessing this defaults to 'densitymx' even for a unitarily parameterized model?) we don't encounter the error. This also doesn't look to be related to the use of a statically parameterized model, as the above behavior is reproduced when using 'full unitary' as the parameterization as well (but again only when specifying the 'statevec' evotype).  I've also inspected the state prep model members directly and can confirm that the '_basis' attribute is None in the cases where this fails.

**Environment:**
 - pyGSTi version 0.9.11.2
 - python version 3.9.13
 - OS: Windows 11



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index e89bbe20f..b2cbd34b7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -47,6 +47,8 @@ jupyter_notebooks/Tutorials/tutorial_files/exampleMultiDataSetReport
 jupyter_notebooks/Tutorials/tutorial_files/exampleBriefReport
 jupyter_notebooks/Tutorials/tutorial_files/*.ipynb
 jupyter_notebooks/Tutorials/tutorial_files/tempTest
+jupyter_notebooks/Tutorials/tutorial_files/*checkpoints
+
 
 
 # Compiled source #
diff --git a/jupyter_notebooks/Examples/BootstrappedErrorBars.ipynb b/jupyter_notebooks/Examples/BootstrappedErrorBars.ipynb
index 5e5c77f94..ddf48e4b4 100644
--- a/jupyter_notebooks/Examples/BootstrappedErrorBars.ipynb
+++ b/jupyter_notebooks/Examples/BootstrappedErrorBars.ipynb
@@ -49,8 +49,8 @@
     "\n",
     "\n",
     "results = pygsti.run_stdpractice_gst(ds, target_model, prep_fiducials, meas_fiducials,\n",
-    "                                    germs, maxLengths, modes=\"TP\")\n",
-    "estimated_model = results.estimates['TP'].models['stdgaugeopt']"
+    "                                    germs, maxLengths, modes=\"full TP\")\n",
+    "estimated_model = results.estimates['full TP'].models['stdgaugeopt']"
    ]
   },
   {
diff --git a/jupyter_notebooks/Examples/CirqIntegration.ipynb b/jupyter_notebooks/Examples/CirqIntegration.ipynb
index 8ebfe8cb4..c5bab68d0 100644
--- a/jupyter_notebooks/Examples/CirqIntegration.ipynb
+++ b/jupyter_notebooks/Examples/CirqIntegration.ipynb
@@ -478,7 +478,7 @@
     }
    ],
    "source": [
-    "gst_results = pygsti.run_stdpractice_gst(dataset, target_model, preps, effects, germs, max_lengths, modes=\"TP,Target\", verbosity=1)"
+    "gst_results = pygsti.run_stdpractice_gst(dataset, target_model, preps, effects, germs, max_lengths, modes=[\"full TP\",\"Target\"], verbosity=1)"
    ]
   },
   {
@@ -503,7 +503,7 @@
     }
    ],
    "source": [
-    "mdl_estimate = gst_results.estimates['TP'].models['stdgaugeopt']\n",
+    "mdl_estimate = gst_results.estimates['full TP'].models['stdgaugeopt']\n",
     "print(\"2DeltaLogL(estimate, data): \", pygsti.tools.two_delta_logl(mdl_estimate, dataset))\n",
     "print(\"2DeltaLogL(ideal, data): \", pygsti.tools.two_delta_logl(target_model, dataset))"
    ]
diff --git a/jupyter_notebooks/Examples/GOpt-AddingNewOptimizations.ipynb b/jupyter_notebooks/Examples/GOpt-AddingNewOptimizations.ipynb
index 7df2cae38..b02e93753 100644
--- a/jupyter_notebooks/Examples/GOpt-AddingNewOptimizations.ipynb
+++ b/jupyter_notebooks/Examples/GOpt-AddingNewOptimizations.ipynb
@@ -32,7 +32,7 @@
     "ds = pygsti.data.simulate_data(mdl_datagen, exp_design.all_circuits_needing_data, num_samples=1000, seed=1234)\n",
     "data = pygsti.protocols.ProtocolData(exp_design, ds)\n",
     "\n",
-    "gst = pygsti.protocols.StandardGST(\"TP\", gaugeopt_suite={'go0': {'item_weights': {'gates': 1, 'spam': 1}}})\n",
+    "gst = pygsti.protocols.StandardGST(\"full TP\", gaugeopt_suite={'go0': {'item_weights': {'gates': 1, 'spam': 1}}})\n",
     "results = gst.run(data) \n",
     "results.write(\"example_files/regaugeopt_example\")"
    ]
@@ -59,7 +59,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "estimate = my_results.estimates['TP']\n",
+    "estimate = my_results.estimates['full TP']\n",
     "estimate.add_gaugeoptimized( {'item_weights': {'gates': 1, 'spam': 0.001}}, label=\"Spam 1e-3\" )\n",
     "mdl_gaugeopt = estimate.models['Spam 1e-3']\n",
     "\n",
diff --git a/jupyter_notebooks/Examples/GOpt-NonIdealTargets.ipynb b/jupyter_notebooks/Examples/GOpt-NonIdealTargets.ipynb
index 874773977..da4d84172 100644
--- a/jupyter_notebooks/Examples/GOpt-NonIdealTargets.ipynb
+++ b/jupyter_notebooks/Examples/GOpt-NonIdealTargets.ipynb
@@ -72,7 +72,7 @@
    "outputs": [],
    "source": [
     "# GST with standard \"ideal target\" gauge optimization\n",
-    "results1 = pygsti.protocols.StandardGST(\"TP\").run(data)"
+    "results1 = pygsti.protocols.StandardGST(\"full TP\").run(data)"
    ]
   },
   {
@@ -84,7 +84,7 @@
     "# GST with our guess as the gauge optimization target\n",
     "gaugeopt_suite = pygsti.protocols.GSTGaugeOptSuite(gaugeopt_suite_names=['stdgaugeopt'],\n",
     "                                                  gaugeopt_target=mdl_guess)\n",
-    "results2 = pygsti.protocols.StandardGST(\"TP\", gaugeopt_suite).run(data)"
+    "results2 = pygsti.protocols.StandardGST(\"full TP\", gaugeopt_suite).run(data)"
    ]
   },
   {
@@ -104,8 +104,8 @@
    "outputs": [],
    "source": [
     "target_model = smq1Q_XYI.target_model()\n",
-    "mdl_1 = results1.estimates['TP'].models['stdgaugeopt']\n",
-    "mdl_2 = results2.estimates['TP'].models['stdgaugeopt']\n",
+    "mdl_1 = results1.estimates['full TP'].models['stdgaugeopt']\n",
+    "mdl_2 = results2.estimates['full TP'].models['stdgaugeopt']\n",
     "print(\"Diff between ideal and ideal-target-gauge-opt = \", mdl_1.frobeniusdist(target_model))\n",
     "print(\"Diff between ideal and mdl_guess-gauge-opt = \", mdl_2.frobeniusdist(target_model))\n",
     "print(\"Diff between data-gen and ideal-target-gauge-opt = \", mdl_1.frobeniusdist(mdl_datagen))\n",
@@ -133,7 +133,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results1.estimates['TP'].add_gaugeoptimized(results2.estimates['TP'].goparameters['stdgaugeopt'],\n",
+    "results1.estimates['full TP'].add_gaugeoptimized(results2.estimates['full TP'].goparameters['stdgaugeopt'],\n",
     "                                            label=\"using mdl_guess\")"
    ]
   },
@@ -143,7 +143,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "mdl_1b = results1.estimates['TP'].models['using mdl_guess']\n",
+    "mdl_1b = results1.estimates['full TP'].models['using mdl_guess']\n",
     "print(mdl_1b.frobeniusdist(mdl_2)) # gs1b is the same as gs2"
    ]
   },
diff --git a/jupyter_notebooks/Examples/MPI-RunningGST.ipynb b/jupyter_notebooks/Examples/MPI-RunningGST.ipynb
index 227cce87f..44b2ab235 100644
--- a/jupyter_notebooks/Examples/MPI-RunningGST.ipynb
+++ b/jupyter_notebooks/Examples/MPI-RunningGST.ipynb
@@ -62,7 +62,7 @@
     "memLim = 2.1*(1024)**3  # 2.1 GB\n",
     "\n",
     "#Perform TP-constrained GST\n",
-    "protocol = pygsti.protocols.StandardGST(\"TP\")\n",
+    "protocol = pygsti.protocols.StandardGST(\"full TP\")\n",
     "start = time.time()\n",
     "results = protocol.run(data, memlimit=memLim, comm=comm)\n",
     "end = time.time()\n",
diff --git a/jupyter_notebooks/Tutorials/algorithms/GST-Overview.ipynb b/jupyter_notebooks/Tutorials/algorithms/GST-Overview.ipynb
index 986cb3c8a..932b5fd88 100644
--- a/jupyter_notebooks/Tutorials/algorithms/GST-Overview.ipynb
+++ b/jupyter_notebooks/Tutorials/algorithms/GST-Overview.ipynb
@@ -106,7 +106,7 @@
    "outputs": [],
    "source": [
     "#run the GST protocol and create a report \n",
-    "gst_protocol = pygsti.protocols.StandardGST('full TP,CPTP,Target')\n",
+    "gst_protocol = pygsti.protocols.StandardGST(['full TP','CPTPLND','Target'])\n",
     "results = gst_protocol.run(data)\n",
     "\n",
     "report = pygsti.report.construct_standard_report(\n",
diff --git a/jupyter_notebooks/Tutorials/algorithms/GST-Protocols.ipynb b/jupyter_notebooks/Tutorials/algorithms/GST-Protocols.ipynb
index 8427880c6..0ab42b10a 100644
--- a/jupyter_notebooks/Tutorials/algorithms/GST-Protocols.ipynb
+++ b/jupyter_notebooks/Tutorials/algorithms/GST-Protocols.ipynb
@@ -28,7 +28,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from __future__ import print_function\n",
     "import pygsti"
    ]
   },
@@ -63,7 +62,7 @@
     "## `GateSetTomography`\n",
     "This protocol performs a single model optimization, and so computes a **single GST estimate** given a `DataSet`, a target `Model`, and other parameters.  (The returned `ModelEstimateResults` object may sometimes contain multiple related estimates in certain cases, but in these cases all the estimates are closely related.)  The experiment design provides all of the information about the GST circuits, in this case a *standard*  (*prep_fiducial + germ^power + meas_fiducial*) set, so the only thing needed by the protocol is an initial `Model` to optimize.  Thus, the `GateSetTomography` protocol is essentially just a model optimizer that you give an initial point.  Importantly, this initial point (a `Model`) also specifies the *parameterization*, i.e. the space of parameters that are optimized over.\n",
     "\n",
-    "Minimally, when using `GateSetTomography` you should set the parameterization of the initial model.  This can be viewed as setting the constraints on the optimization.  For instance, when the gates in the model are parameterized as trace-preserving (TP) maps, the optimization will be constrained to trying gate sets with TP gates (because every set of parameters corresponds to a set of TP gates).  In the cell below, we constrain the optimization to TP gate sets by using `.target_model(\"TP\")`, which returns a version of the target model where all the gates are TP-parameterized, the state preparation has trace = 1, and the POVM effects always add to the identity.  This could also be done by calling `set_all_parameterizations(\"TP\")` on the fully-parameterized target model returned by `.target_model()`.  See the [tutorial on explicit models](../objects/ExplicitModel.ipynb) for more information on setting a model's parameterization."
+    "Minimally, when using `GateSetTomography` you should set the parameterization of the initial model.  This can be viewed as setting the constraints on the optimization.  For instance, when the gates in the model are parameterized as trace-preserving (TP) maps, the optimization will be constrained to trying gate sets with TP gates (because every set of parameters corresponds to a set of TP gates).  In the cell below, we constrain the optimization to TP gate sets by using `.target_model(\"full TP\")`, which returns a version of the target model where all the gates are TP-parameterized, the state preparation has trace = 1, and the POVM effects always add to the identity.  This could also be done by calling `set_all_parameterizations(\"TP\")` on the fully-parameterized target model returned by `.target_model()`.  See the [tutorial on explicit models](../objects/ExplicitModel.ipynb) for more information on setting a model's parameterization."
    ]
   },
   {
@@ -126,7 +125,7 @@
     "    target_model_TP2, name=\"GSTwithMyGO\",\n",
     "    gaugeopt_suite={'my_gauge_opt': {'item_weights': {'gates': 1.0, 'spam': 0.001}}}\n",
     "    )\n",
-    "results_TP2 = proto.run(data)"
+    "results_TP2 = proto.run(data, disable_checkpointing=True)"
    ]
   },
   {
@@ -168,7 +167,7 @@
     "\n",
     "\n",
     "proto = pygsti.protocols.GateSetTomography(target_model_TP2, name=\"GSTwithReducedData\")\n",
-    "results_reduced = proto.run(reduced_data)"
+    "results_reduced = proto.run(reduced_data, disable_checkpointing=True)"
    ]
   },
   {
@@ -183,7 +182,7 @@
    "metadata": {},
    "source": [
     "## `StandardGST`\n",
-    "The protocol embodies a standard *set* of GST protocols to be run on a set of data.  It essentially runs multiple `GateSetTomography` protocols on the given data which use different parameterizations of an `ExplicitOpModel`  (the `StandardGST` protocol doesn't work with other types of `Model` objects, e.g. *implicit* models, which don't implement `set_all_parameterizations`).  The `modes` argument is a comma-separated list of the parameterization types that should be run (e.g. `\"TP,CPTP\"` will compute a Trace-Preserving estimate *and* a Completely-Positive & Trace-Preserving estimate). The currently available modes are:\n",
+    "The protocol embodies a standard *set* of GST protocols to be run on a set of data.  It essentially runs multiple `GateSetTomography` protocols on the given data which use different parameterizations of an `ExplicitOpModel`  (the `StandardGST` protocol doesn't work with other types of `Model` objects, e.g. *implicit* models, which don't implement `set_all_parameterizations`).  The `modes` argument is a list strings corresponding to the parameterization types that should be run (e.g. `[\"full TP\",\"CPTPLND\"]` will compute a Trace-Preserving estimate *and* a Completely-Positive & Trace-Preserving estimate). The currently available modes are:\n",
     " - \"full\" : unconstrained gates (fully parameterized)                                                                 \n",
     " - \"TP\"   : TP-constrained gates and state preparations\n",
     " - \"CPTP\" : CPTP-constrained gates and TP-constrained state preparations               \n",
@@ -321,6 +320,163 @@
     "# pickle.dump(results_stdprac, open('../tutorial_files/exampleResults_stdprac.pkl',\"wb\"))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Checkpointing/Warmstarting\n",
+    "\n",
+    "The `GateSetTomography` and `StandardGST` protocols both support checkpointing to enable resuming GST analysis after an unexpected failure, such as an out-of-memory error, or an unexpected timeout in resource limited compute environments (clusters etc.), or for whatever other reason. Checkpointing is enabled by default, so no additional changes are needed in order to have these generated. \n",
+    "\n",
+    "Each protocol has a corresponding checkpoint object, `GateSetTomographyCheckpoint` and `StandardGSTCheckpoint`, which are saved to disk over the course of an iterative fit in serialized json format. By default checkpoint files associated with a `GateSetTomographyCheckpoint` object are saved to a new directory located in whichever current working directory the protocol is being run from named 'gst_checkpoints'. A new file is written to disk after each iteration with default naming of the form `GateSetTomography_iteration_{i}.json` where i is the index of the completed GST iteration associated with that checkpoint. Similarly, for a `StandardGSTCheckpoint` object the checkpoints are by default saved to a directory named 'standard_gst_checkpoints' with default file names of the form `StandardGST_{mode}_iteration_{i}` where mode corresponds to the current parameterized fit or model test associated with that file (including checkpoint information for all previously completed modes prior to the currently running one) and i is the index of the completed iteration within that current mode.\n",
+    "\n",
+    "Below we repeat our first example of the notebook, but this time with checkpointing enabled (as is the default)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pygsti.modelpacks import smq1Q_XYI\n",
+    "target_model_TP = smq1Q_XYI.target_model(\"full TP\")\n",
+    "proto = pygsti.protocols.GateSetTomography(target_model_TP)\n",
+    "results_TP = proto.run(data, checkpoint_path = '../tutorial_files/gst_checkpoints/GateSetTomography')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that in the example above we have specified a value for an additional kwarg called `checkpoint_path`. This allows for overriding the default behavior for the save location and naming of checkpoint files. The expected format is `{path}/{name}` where path is the directory to save the checkpoint files to (with that directory being created is required) and where name is the stem of the checkpoint file names `{name}_iteration_{i}.json`. Inspecting the contents of the directory we just specified, we can see that it is now populated by 8 new checkpoint files."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "os.listdir('../tutorial_files/gst_checkpoints/')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Suppose hypothetically that a GST fit had failed at iteration 5 and we wanted to restart from that point without redoing all of the previous iterations from scratch again. We'll call this warmstarting. We can do so by reading in the appropriate serialized checkpoint object using the `read` class method of `GateSetTomographyCheckpoint` and passing that now loaded checkpoint object in for the `checkpoint` kwarg of `run`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pygsti.protocols import GateSetTomographyCheckpoint\n",
+    "gst_iter_5_checkpoint = GateSetTomographyCheckpoint.read('../tutorial_files/gst_checkpoints/GateSetTomography_iteration_5.json')\n",
+    "results_TP_from_iter_5= proto.run(data, checkpoint= gst_iter_5_checkpoint, checkpoint_path = '../tutorial_files/gst_checkpoints/GateSetTomography')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can see from the output that we indeed started from iteration 6 (note the output log indexes from 1 instead of 0). Moreover we can see that we've indeed produced the same output as before without warmstarting, as we would expect/hope:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "all(results_TP.estimates['GateSetTomography'].models['final iteration estimate'].to_vector() == \\\n",
+    "results_TP_from_iter_5.estimates['GateSetTomography'].models['final iteration estimate'].to_vector())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The checkpoint object itself contains information that could be useful for diagnostics or debugging, including the current list of models associated each iterative fit, the last completed iteration it is associated with, and the list of circuits for the last completed iteration it is associated with."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Checkpointing with the StandardGST protocol works similarly:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "proto_standard_gst = pygsti.protocols.StandardGST(modes=['full TP', 'CPTPLND', 'Target'], verbosity=3)\n",
+    "results_stdprac = proto_standard_gst.run(data, checkpoint_path = '../tutorial_files/standard_gst_checkpoints/StandardGST')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Except this time we have significantly more files saved, as during the course of the StandardGST protocol we're actually running three subprotocols:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "os.listdir('../tutorial_files/standard_gst_checkpoints/')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that the StandardGST protocol runs the subprotocols in the order listed in the `modes` argument, and checkpoint objects labeled with a given model label additionally contain the checkpointing information for the final iterations of any preceding modes which have been completed. i.e. the CPTPLND checkpoint objects contain the information required for full TP. Likewise, checkpoints for Target contain the information required for the full TP and CPTPLND modes. As before, imagine that our fitting failed for whatever reason during iteration 5 of CPTPLND, we can warmstart the protocol by loading in the checkpoint object associated with iteration 4 as below:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "from pygsti.protocols import StandardGSTCheckpoint\n",
+    "standard_gst_checkpoint = StandardGSTCheckpoint.read('../tutorial_files/standard_gst_checkpoints/StandardGST_CPTPLND_iteration_4.json')\n",
+    "results_stdprac_warmstart= proto_standard_gst.run(data, checkpoint= standard_gst_checkpoint, checkpoint_path = '../tutorial_files/standard_gst_checkpoints/StandardGST')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notice that we've indeed skipped past the previously completed full TP mode and jumped straight to the 6th iteration of the CPTPLND fit as expected. \n",
+    "\n",
+    "As for the GateSetTomographyCheckpoint object described above, the `StandardGSTCheckpoint` can often be useful to inspect as a debugging/diagnostic tool. `StandardGSTCheckpoints` are essentially structured as container object that hold a set of child `GateSetTomographyCheckpoint` and `ModelTestCheckpoint` (more on these in the ModelTest tutorial) objects for each of the modes being run (and potentially more types of chile checkpoints in the future as we add additional functionality). These children can be accessed using the `children` attribute of a `StandardGSTCheckpoint` instance which is a dictionary with keys given by the mode names contained therein."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(standard_gst_checkpoint.children['CPTPLND'])"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -331,9 +487,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "gst_checkpointing",
    "language": "python",
-   "name": "python3"
+   "name": "gst_checkpointing"
   },
   "language_info": {
    "codemirror_mode": {
diff --git a/jupyter_notebooks/Tutorials/algorithms/ModelTesting.ipynb b/jupyter_notebooks/Tutorials/algorithms/ModelTesting.ipynb
index b8166a8b0..1fe33abc5 100644
--- a/jupyter_notebooks/Tutorials/algorithms/ModelTesting.ipynb
+++ b/jupyter_notebooks/Tutorials/algorithms/ModelTesting.ipynb
@@ -13,12 +13,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "from __future__ import division, print_function\n",
-    "\n",
     "import pygsti\n",
     "import numpy as np\n",
     "import scipy\n",
@@ -28,7 +26,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -49,7 +47,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -72,100 +70,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# creates a Results object with a \"model1\" estimate\n",
     "results = pygsti.protocols.ModelTest(test_model1, target_model, name='model1').run(data)\n",
@@ -186,21 +93,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Running idle tomography\n",
-      "Computing switchable properties\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "results.add_estimates(results2)\n",
     "results.add_estimates(results3)\n",
@@ -229,57 +124,9 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "-- Std Practice:  Iter 1 of 3  (full TP) --: \n",
-      "  --- Iterative GST: [##################################################] 100.0%  952 circuits ---\n",
-      "  Iterative GST Total Time: 3.0s\n",
-      "-- Std Practice:  Iter 2 of 3  (CPTP) --: \n",
-      "  --- Iterative GST: [##################################################] 100.0%  952 circuits ---\n",
-      "  Iterative GST Total Time: 6.7s\n",
-      "-- Std Practice:  Iter 3 of 3  (Target) --: \n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "  MatrixLayout: 1 processors divided into 1 (= 1) grid along circuit and parameter directions.\n",
-      "     1 atoms, parameter block size limits ()\n",
-      "  *** Distributing 1 atoms to 1 atom-processing groups (1 cores) ***\n",
-      "      More atom-processors than hosts: each host gets ~1 atom-processors\n",
-      "Running idle tomography\n",
-      "Computing switchable properties\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "#Create some GST results using standard practice GST\n",
     "gst_results = pygsti.protocols.StandardGST().run(data)\n",
@@ -307,32 +154,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "-- Std Practice:  Iter 1 of 4  (full TP) --: \n",
-      "  --- Iterative GST: [##################################################] 100.0%  952 circuits ---\n",
-      "  Iterative GST Total Time: 3.2s\n",
-      "-- Std Practice:  Iter 2 of 4  (Test2) --: \n",
-      "-- Std Practice:  Iter 3 of 4  (Test3) --: \n",
-      "-- Std Practice:  Iter 4 of 4  (Target) --: \n",
-      "Running idle tomography\n",
-      "Computing switchable properties\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n",
-      "Found standard clifford compilation from smq1Q_XYI\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "proto = pygsti.protocols.StandardGST(modes=\"full TP,Test2,Test3,Target\", # You MUST put Test2 and Test3 here...\n",
+    "proto = pygsti.protocols.StandardGST(modes=[\"full TP\",\"Test2\",\"Test3\",\"Target\"], # You MUST put Test2 and Test3 here...\n",
     "                                     models_to_test={'Test2': test_model2, 'Test3': test_model3})\n",
-    "gst_results = proto.run(data)\n",
+    "gst_results = proto.run(data, disable_checkpointing=True)\n",
     "\n",
     "pygsti.report.construct_standard_report(\n",
     "    gst_results, title=\"GST with Model Test Example Report 2\", verbosity=1\n",
@@ -343,7 +171,62 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Thats it!  Now that you know more about model-testing you may want to go back to the [overview of pyGST protocls](../03-Protocols.ipynb)."
+    "## Checkpointing/Warmstarting\n",
+    "\n",
+    "Just like the GST protocols discussed in [GST-Protocols](GST-Protocols.ipynb), `ModelTest` protocols also support checkpointing/warmstarting in the event of unexpected failure or termination using `ModelTestCheckpoint` objects which are serialized and written to disk over the course of running an iterative `ModelTest`. We direct you to the linked tutorial for more on the basic syntax and usage, which essentially identical for `ModelTestCheckpoints` objects."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# creates a Results object with a \"model1\" estimate\n",
+    "modeltest_proto = pygsti.protocols.ModelTest(test_model1, target_model, name='model1')\n",
+    "results = modeltest_proto.run(data, checkpoint_path = '../tutorial_files/modeltest_checkpoints/ModelTest_model1')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We now have a bunch of checkpoint files associated with this `ModelTest` protocol saved in the specified directory:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "os.listdir('../tutorial_files/modeltest_checkpoints/')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Suppose hypothetically we needed to restart the running of this ModelTest from where it left off following iteration 4. We can do so reading in the saved `ModelTestCheckpoint` object and passing it into the `ModelTest`'s `run` method."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pygsti.protocols import ModelTestCheckpoint\n",
+    "modeltest_checkpoint = ModelTestCheckpoint.read('../tutorial_files/modeltest_checkpoints/ModelTest_model1_iteration_4.json')\n",
+    "results_warmstarted= modeltest_proto.run(data, checkpoint= modeltest_checkpoint, checkpoint_path= '../tutorial_files/modeltest_checkpoints/ModelTest_model1')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As expected, the results are identical in both cases:"
    ]
   },
   {
@@ -351,14 +234,25 @@
    "execution_count": null,
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "print(results.estimates['model1'].parameters['model_test_values'])\n",
+    "print(results.estimates['model1'].parameters['model_test_values'] == \\\n",
+    "results_warmstarted.estimates['model1'].parameters['model_test_values'])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Thats it!  Now that you know more about model-testing you may want to go back to the [overview of pyGSTi protocols](../03-Protocols.ipynb)."
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "gst_checkpointing",
    "language": "python",
-   "name": "python3"
+   "name": "gst_checkpointing"
   },
   "language_info": {
    "codemirror_mode": {
diff --git a/jupyter_notebooks/Tutorials/reporting/ReportGeneration.ipynb b/jupyter_notebooks/Tutorials/reporting/ReportGeneration.ipynb
index e7c8ab329..4b8facc59 100644
--- a/jupyter_notebooks/Tutorials/reporting/ReportGeneration.ipynb
+++ b/jupyter_notebooks/Tutorials/reporting/ReportGeneration.ipynb
@@ -151,7 +151,7 @@
    "source": [
     "ws = pygsti.report.Workspace()\n",
     "report = pygsti.report.construct_standard_report(\n",
-    "    {'TP': results_tp, \"Full\": results_full}, title=\"Example Multi-Estimate Report\", ws=ws, verbosity=2)\n",
+    "    {'full TP': results_tp, \"Full\": results_full}, title=\"Example Multi-Estimate Report\", ws=ws, verbosity=2)\n",
     "report.write_html(\"../tutorial_files/exampleMultiEstimateReport\", auto_open=True, verbosity=2)"
    ]
   },
diff --git a/jupyter_notebooks/Tutorials/reporting/WorkspaceExamples.ipynb b/jupyter_notebooks/Tutorials/reporting/WorkspaceExamples.ipynb
index f0a5ad7eb..2d4b3c2c7 100644
--- a/jupyter_notebooks/Tutorials/reporting/WorkspaceExamples.ipynb
+++ b/jupyter_notebooks/Tutorials/reporting/WorkspaceExamples.ipynb
@@ -265,7 +265,7 @@
     "             [4,(1.0,0.5,0)],[10,(1.0,0,0)]])\n",
     "#print(cMap.colorscale())\n",
     "w.MatrixPlot(mx, colormap=cMap, colorbar=False, grid=\"white:1\", box_labels=True, prec=2,\n",
-    "             xlabels=('TP',\"CPTP\",\"full\"),ylabels=(\"DS0\",\"DS1\",\"DS2\"))"
+    "             xlabels=('full TP',\"CPTPLND\",\"full\"),ylabels=(\"DS0\",\"DS1\",\"DS2\"))"
    ]
   },
   {
diff --git a/pygsti/algorithms/core.py b/pygsti/algorithms/core.py
index f44125ca9..640b4174c 100644
--- a/pygsti/algorithms/core.py
+++ b/pygsti/algorithms/core.py
@@ -749,6 +749,89 @@ def run_iterative_gst(dataset, start_model, circuit_lists,
     final_objfn : MDSObjectiveFunction
         The final iteration's objective function / store, which encapsulated the final objective
         function evaluated at the best-fit point (an "evaluated" model-dataSet-circuits store).
+
+    """
+    gst_iter_gen = iterative_gst_generator(dataset, start_model, circuit_lists,
+                      optimizer, iteration_objfn_builders, final_objfn_builders,
+                      resource_alloc, starting_index=0, verbosity=verbosity)
+
+    models = []
+    optimums = []
+
+    for i in range(len(circuit_lists)):
+        #then do the final iteration slightly differently since the generator should
+        #give three return values.
+        if i==len(circuit_lists)-1:
+            mdl_iter, opt_iter, final_objfn =  next(gst_iter_gen)
+        else:
+            mdl_iter, opt_iter =  next(gst_iter_gen)
+
+        models.append(mdl_iter)
+        optimums.append(opt_iter)
+        
+    return models, optimums, final_objfn
+
+def iterative_gst_generator(dataset, start_model, circuit_lists,
+                      optimizer, iteration_objfn_builders, final_objfn_builders,
+                      resource_alloc, starting_index=0, verbosity=0):
+    """
+    Performs Iterative Gate Set Tomography on the dataset.
+    Same as `run_iterative_gst`, except this function produces a
+    generator for producing the output for each iteration instead
+    of returning the lists of outputs all at once.
+
+    Parameters
+    ----------
+    dataset : DataSet
+        The data used to generate MLGST gate estimates
+
+    start_model : Model
+        The Model used as a starting point for the least-squares
+        optimization.
+
+    circuit_lists : list of lists of (tuples or Circuits)
+        The i-th element is a list of the circuits to be used in the i-th iteration
+        of the optimization.  Each element of these lists is a circuit, specifed as
+        either a Circuit object or as a tuple of operation labels (but all must be specified
+        using the same type).
+        e.g. [ [ (), ('Gx',) ], [ (), ('Gx',), ('Gy',) ], [ (), ('Gx',), ('Gy',), ('Gx','Gy') ]  ]
+
+    optimizer : Optimizer or dict
+        The optimizer to use, or a dictionary of optimizer parameters
+        from which a default optimizer can be built.
+
+    iteration_objfn_builders : list
+        List of ObjectiveFunctionBuilder objects defining which objective functions
+        should be optimizized (successively) on each iteration.
+
+    final_objfn_builders : list
+        List of ObjectiveFunctionBuilder objects defining which objective functions
+        should be optimizized (successively) on the final iteration.
+    
+    resource_alloc : ResourceAllocation
+        A resource allocation object containing information about how to
+        divide computation amongst multiple processors and any memory
+        limits that should be imposed.
+
+    
+    starting_index : int, optional (default 0)
+        Index of the iteration to start the optimization at. Primarily used
+        when warmstarting the iterative optimization from a checkpoint.
+    
+    verbosity : int, optional
+        How much detail to send to stdout.
+
+
+    Returns
+    -------
+    generator
+        Returns a generator which when queried the i-th time returns a tuple containing: 
+
+        - model: the model corresponding to the results of the i-th iteration.
+        - optimums : the final OptimizerResults from the i-th iteration.
+        - final_objfn : If the final iteration the MDSObjectiveFunction function / store, 
+          which encapsulated the final objective function evaluated at the best-fit point 
+          (an "evaluated" model-dataset-circuits store).
     """
     resource_alloc = _ResourceAllocation.cast(resource_alloc)
     optimizer = optimizer if isinstance(optimizer, _Optimizer) else _CustomLMOptimizer.cast(optimizer)
@@ -756,12 +839,10 @@ def run_iterative_gst(dataset, start_model, circuit_lists,
     profiler = resource_alloc.profiler
     printer = VerbosityPrinter.create_printer(verbosity, comm)
 
-    models = []; optimums = []
     mdl = start_model.copy(); nIters = len(circuit_lists)
     tStart = _time.time()
     tRef = tStart
     final_objfn = None
-
     iteration_objfn_builders = [_objfns.ObjectiveFunctionBuilder.cast(ofb) for ofb in iteration_objfn_builders]
     final_objfn_builders = [_objfns.ObjectiveFunctionBuilder.cast(ofb) for ofb in final_objfn_builders]
 
@@ -777,8 +858,10 @@ def _max_array_types(artypes_list):  # get the maximum number of each array type
         return ret
 
     with printer.progress_logging(1):
-        for (i, circuitsToEstimate) in enumerate(circuit_lists):
+        for i in range(starting_index, len(circuit_lists)):
+            circuitsToEstimate = circuit_lists[i]
             extraMessages = []
+
             if isinstance(circuitsToEstimate, _CircuitList) and circuitsToEstimate.name:
                 extraMessages.append("(%s) " % circuitsToEstimate.name)
 
@@ -802,6 +885,7 @@ def _max_array_types(artypes_list):  # get the maximum number of each array type
                     first_iter_optimizer = _copy.deepcopy(optimizer)  # use a separate copy of optimizer, as it
                     first_iter_optimizer.fditer = optimizer.first_fditer  # is a persistent object (so don't modify!)
                     opt_result, mdc_store = run_gst_fit(mdc_store, first_iter_optimizer, obj_fn_builder, printer - 1)
+
                 else:
                     opt_result, mdc_store = run_gst_fit(mdc_store, optimizer, obj_fn_builder, printer - 1)
                 profiler.add_time('run_iterative_gst: iter %d %s-opt' % (i + 1, obj_fn_builder.name), tNxt)
@@ -823,20 +907,19 @@ def _max_array_types(artypes_list):  # get the maximum number of each array type
                 printer.log("Final optimization took %.1fs\n" % (tNxt - tRef), 2)
                 tRef = tNxt
 
-                models.append(mdc_store.model)  # don't copy so `mdc_store.model` *is* the final model, `models[-1]`
-
+                # don't copy so `mdc_store.model` *is* the final model, `models[-1]`
                 # send final objfn object back to caller to facilitate postproc  on the final (model, circuits, dataset)
                 # Note: initial_mdc_store is *not* an objective fn (it's just a store) so don't send it back.
                 if mdc_store is not initial_mdc_store:
                     final_objfn = mdc_store
-            else:
-                models.append(mdc_store.model.copy())
 
-            optimums.append(opt_result)
+                yield (mdc_store.model, opt_result, final_objfn)
+            else:
+                #If not the final iteration then only send back a copy of the model and the optimizer results
+                yield (mdc_store.model.copy(), opt_result)
 
     printer.log('Iterative GST Total Time: %.1fs' % (_time.time() - tStart))
     profiler.add_time('run_iterative_gst: total time', tStart)
-    return models, optimums, final_objfn
 
 
 def _do_runopt(objective, optimizer, printer):
diff --git a/pygsti/drivers/longsequence.py b/pygsti/drivers/longsequence.py
index 423d08061..b1caabd0b 100644
--- a/pygsti/drivers/longsequence.py
+++ b/pygsti/drivers/longsequence.py
@@ -34,7 +34,8 @@ def run_model_test(model_filename_or_object,
                    prep_fiducial_list_or_filename, meas_fiducial_list_or_filename,
                    germs_list_or_filename, max_lengths, gauge_opt_params=None,
                    advanced_options=None, comm=None, mem_limit=None,
-                   output_pkl=None, verbosity=2):
+                   output_pkl=None, verbosity=2, checkpoint=None, checkpoint_path=None,
+                   disable_checkpointing= False):
     """
     Compares a :class:`Model`'s predictions to a `DataSet` using GST-like circuits.
 
@@ -120,6 +121,23 @@ def run_model_test(model_filename_or_object,
         The 'verbosity' option is an integer specifying the level of
         detail printed to stdout during the calculation.
 
+    checkpoint : ModelTestCheckpoint, optional (default None)
+            If specified use a previously generated checkpoint object to restart
+            or warm start this run part way through.
+
+    checkpoint_path : str, optional (default None)
+        A string for the path/name to use for writing intermediate checkpoint
+        files to disk. Format is {path}/{name}, without inclusion of the json
+        file extension. This {path}/{name} combination will have the latest
+        completed iteration number appended to it before writing it to disk.
+        If none, the value of {name} will be set to the name of the protocol
+        being run.
+
+    disable_checkpointing : bool, optional (default False)
+        When set to True checkpoint objects will not be constructed and written
+        to disk during the course of this protocol. It is strongly recommended
+        that this be kept set to False without good reason to disable the checkpoints.
+
     Returns
     -------
     Results
@@ -167,7 +185,7 @@ def run_model_test(model_filename_or_object,
     proto.circuit_weights = advanced_options.get('circuit_weights', None)
     proto.unreliable_ops = advanced_options.get('unreliable_ops', ['Gcnot', 'Gcphase', 'Gms', 'Gcn', 'Gcx', 'Gcz'])
 
-    results = proto.run(data, mem_limit, comm)
+    results = proto.run(data, mem_limit, comm, checkpoint=checkpoint, checkpoint_path=checkpoint_path, disable_checkpointing=disable_checkpointing)
     _output_to_pickle(results, output_pkl, comm)
     return results
 
@@ -287,7 +305,8 @@ def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
                           prep_fiducial_list_or_filename, meas_fiducial_list_or_filename,
                           germs_list_or_filename, max_lengths, gauge_opt_params=None,
                           advanced_options=None, comm=None, mem_limit=None,
-                          output_pkl=None, verbosity=2):
+                          output_pkl=None, verbosity=2, checkpoint=None, checkpoint_path=None,
+                          disable_checkpointing = False):
     """
     Perform long-sequence GST (LSGST).
 
@@ -409,6 +428,22 @@ def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
         - 4 -- also shows inner iterations of LM algorithm
         - 5 -- also shows detailed info from within jacobian and objective function calls.
 
+    checkpoint : GateSetTomographyCheckpoint, optional (default None)
+        If specified use a previously generated checkpoint object to restart
+        or warm start this run part way through.
+
+    checkpoint_path : str, optional (default None)
+        A string for the path/name to use for writing intermediate checkpoint
+        files to disk. Format is {path}/{name}, without inclusion of the json
+        file extension. This {path}/{name} combination will have the latest
+        completed iteration number appended to it before writing it to disk.
+        If none, the value of {name} will be set to the name of the protocol
+        being run.
+    disable_checkpointing : bool, optional (default False)
+        When set to True checkpoint objects will not be constructed and written
+        to disk during the course of this protocol. It is strongly recommended
+        that this be kept set to False without good reason to disable the checkpoints.
+
     Returns
     -------
     Results
@@ -453,7 +488,7 @@ def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
     proto.circuit_weights = advanced_options.get('circuit_weights', None)
     proto.unreliable_ops = advanced_options.get('unreliable_ops', ['Gcnot', 'Gcphase', 'Gms', 'Gcn', 'Gcx', 'Gcz'])
 
-    results = proto.run(data, mem_limit, comm)
+    results = proto.run(data, mem_limit, comm, checkpoint=checkpoint, checkpoint_path= checkpoint_path, disable_checkpointing=disable_checkpointing)
     _output_to_pickle(results, output_pkl, comm)
     return results
 
@@ -461,7 +496,8 @@ def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
 def run_long_sequence_gst_base(data_filename_or_set, target_model_filename_or_object,
                                lsgst_lists, gauge_opt_params=None,
                                advanced_options=None, comm=None, mem_limit=None,
-                               output_pkl=None, verbosity=2):
+                               output_pkl=None, verbosity=2, checkpoint=None, checkpoint_path=None,
+                               disable_checkpointing = False):
     """
     A more fundamental interface for performing end-to-end GST.
 
@@ -530,6 +566,23 @@ def run_long_sequence_gst_base(data_filename_or_set, target_model_filename_or_ob
         - 4 -- also shows inner iterations of LM algorithm
         - 5 -- also shows detailed info from within jacobian and objective function calls.
 
+    checkpoint : GateSetTomographyCheckpoint, optional (default None)
+        If specified use a previously generated checkpoint object to restart
+        or warm start this run part way through.
+
+    checkpoint_path : str, optional (default None)
+        A string for the path/name to use for writing intermediate checkpoint
+        files to disk. Format is {path}/{name}, without inclusion of the json
+        file extension. This {path}/{name} combination will have the latest
+        completed iteration number appended to it before writing it to disk.
+        If none, the value of {name} will be set to the name of the protocol
+        being run.
+
+    disable_checkpointing : bool, optional (default False)
+        When set to True checkpoint objects will not be constructed and written
+        to disk during the course of this protocol. It is strongly recommended
+        that this be kept set to False without good reason to disable the checkpoints.
+        
     Returns
     -------
     Results
@@ -562,16 +615,16 @@ def run_long_sequence_gst_base(data_filename_or_set, target_model_filename_or_ob
     proto.circuit_weights = advanced_options.get('circuit_weights', None)
     proto.unreliable_ops = advanced_options.get('unreliable_ops', ['Gcnot', 'Gcphase', 'Gms', 'Gcn', 'Gcx', 'Gcz'])
 
-    results = proto.run(data, mem_limit, comm)
+    results = proto.run(data, mem_limit, comm, checkpoint=checkpoint, checkpoint_path=checkpoint_path, disable_checkpointing=disable_checkpointing)
     _output_to_pickle(results, output_pkl, comm)
     return results
 
 
 def run_stdpractice_gst(data_filename_or_set, target_model_filename_or_object, prep_fiducial_list_or_filename,
                         meas_fiducial_list_or_filename, germs_list_or_filename, max_lengths,
-                        modes="full TP,CPTP,Target", gaugeopt_suite='stdgaugeopt', gaugeopt_target=None,
+                        modes=('full TP','CPTPLND','Target'), gaugeopt_suite='stdgaugeopt', gaugeopt_target=None,
                         models_to_test=None, comm=None, mem_limit=None, advanced_options=None, output_pkl=None,
-                        verbosity=2):
+                        verbosity=2, checkpoint=None, checkpoint_path=None, disable_checkpointing = False):
     """
     Perform end-to-end GST analysis using standard practices.
 
@@ -612,8 +665,8 @@ def run_stdpractice_gst(data_filename_or_set, target_model_filename_or_object, p
         iteration includes the repeated germs truncated to the L-values *up to*
         and including the i-th one.
 
-    modes : str, optional
-        A comma-separated list of modes which dictate what types of analyses
+    modes : iterable of strs, optional (default ('full TP','CPTPLND','Target')
+        An iterable strings corresponding to modes which dictate what types of analyses
         are performed.  Currently, these correspond to different types of
         parameterizations/constraints to apply to the estimated model.
         The default value is usually fine.  Allowed values are:
@@ -678,6 +731,23 @@ def run_stdpractice_gst(data_filename_or_set, target_model_filename_or_object, p
         The 'verbosity' option is an integer specifying the level of
         detail printed to stdout during the calculation.
 
+    checkpoint : StandardGSTCheckpoint, optional (default None)
+        If specified use a previously generated checkpoint object to restart
+        or warm start this run part way through.
+
+    checkpoint_path : str, optional (default None)
+        A string for the path/name to use for writing intermediate checkpoint
+        files to disk. Format is {path}/{name}, without inclusion of the json
+        file extension. This {path}/{name} combination will have the latest
+        completed iteration number appended to it before writing it to disk.
+        If none, the value of {name} will be set to the name of the protocol
+        being run.
+
+    disable_checkpointing : bool, optional (default False)
+        When set to True checkpoint objects will not be constructed and written
+        to disk during the course of this protocol. It is strongly recommended
+        that this be kept set to False without good reason to disable the checkpoints.
+
     Returns
     -------
     Results
@@ -721,7 +791,7 @@ def run_stdpractice_gst(data_filename_or_set, target_model_filename_or_object, p
                                badfit_options=_get_badfit_options(advanced_options), verbosity=printer,
                                name=advanced_options.get('estimate_label', None))
 
-    results = proto.run(data, mem_limit, comm)
+    results = proto.run(data, mem_limit, comm, checkpoint=checkpoint, checkpoint_path= checkpoint_path, disable_checkpointing=disable_checkpointing)
     _output_to_pickle(results, output_pkl, comm)
     return results
 
diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py
index f65da4fab..d28504149 100644
--- a/pygsti/protocols/gst.py
+++ b/pygsti/protocols/gst.py
@@ -16,6 +16,7 @@
 import pickle as _pickle
 import time as _time
 import warnings as _warnings
+import pathlib as _pathlib
 
 import numpy as _np
 from scipy.stats import chi2 as _chi2
@@ -25,6 +26,7 @@
 from pygsti.protocols.estimate import Estimate as _Estimate
 from pygsti.protocols import protocol as _proto
 from pygsti.protocols.modeltest import ModelTest as _ModelTest
+from pygsti.protocols.modeltest import ModelTestCheckpoint as _ModelTestCheckpoint
 from pygsti import algorithms as _alg
 from pygsti import circuits as _circuits
 from pygsti import io as _io
@@ -41,6 +43,7 @@
 from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
 from pygsti.modelmembers import states as _states, povms as _povms
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
+from pygsti.circuits import Circuit
 
 
 #For results object:
@@ -1265,7 +1268,7 @@ def __init__(self, initial_model=None, gaugeopt_suite='stdgaugeopt',
     #    design = GateSetTomographyDesign(target_model, circuit_lists)
     #    return self.run(_proto.ProtocolData(design, dataset))
 
-    def run(self, data, memlimit=None, comm=None):
+    def run(self, data, memlimit=None, comm=None, checkpoint=None, checkpoint_path=None, disable_checkpointing = False):
         """
         Run this protocol on `data`.
 
@@ -1274,13 +1277,30 @@ def run(self, data, memlimit=None, comm=None):
         data : ProtocolData
             The input data.
 
-        memlimit : int, optional
+        memlimit : int, optional (default None)
             A rough per-processor memory limit in bytes.
 
-        comm : mpi4py.MPI.Comm, optional
+        comm : mpi4py.MPI.Comm, optional (default None)
             When not ``None``, an MPI communicator used to run this protocol
             in parallel.
 
+        checkpoint : GateSetTomographyCheckpoint, optional (default None)
+            If specified use a previously generated checkpoint object to restart
+            or warm start this run part way through.
+
+        checkpoint_path : str, optional (default None)
+            A string for the path/name to use for writing intermediate checkpoint
+            files to disk. Format is {path}/{name}, without inclusion of the json
+            file extension. This {path}/{name} combination will have the latest
+            completed iteration number appended to it before writing it to disk.
+            If none, the value of {name} will be set to the name of the protocol
+            being run.
+
+        disable_checkpointing : bool, optional (default False)
+            When set to True checkpoint objects will not be constructed and written
+            to disk during the course of this protocol. It is strongly recommended
+            that this be kept set to False without good reason to disable the checkpoints.
+
         Returns
         -------
         ModelEstimateResults
@@ -1315,17 +1335,90 @@ def run(self, data, memlimit=None, comm=None):
         tnxt = _time.time(); profiler.add_time('GST: loading', tref); tref = tnxt
         mdl_start = self.initial_model.retrieve_model(data.edesign, self.gaugeopt_suite.gaugeopt_target,
                                                       data.dataset, comm)
+        
+        if not disable_checkpointing:
+            #Set the checkpoint_path variable if None
+            if checkpoint_path is None:
+                checkpoint_path = _pathlib.Path('./gst_checkpoints/' + self.name)
+            else:
+                #cast this to a pathlib path with the file extension (suffix) dropped
+                checkpoint_path = _pathlib.Path(checkpoint_path).with_suffix('')
+            
+            #create the parent directory of the checkpoint if needed:
+            checkpoint_path.parent.mkdir(parents=True, exist_ok=True)
+            
+            #If there is no checkpoint we should start from with the seed model,
+            #otherwise we should seed the next iteration with the last iteration's result.
+            #If there is no checkpoint initialize mdl_lsgst_list and final_objfn to be empty, 
+            #otherwise re-initialize their values from the checkpoint   
+            if checkpoint is None:
+                seed_model = mdl_start.copy()
+                mdl_lsgst_list = []
+                checkpoint = GateSetTomographyCheckpoint()
+            elif isinstance(checkpoint, GateSetTomographyCheckpoint):
+                #if the checkpoint's last completed iteration is non-negative
+                #(i.e. the checkpoint actually has data in it)
+                if checkpoint.last_completed_iter >= 0:
+                    seed_model = checkpoint.mdl_list[-1]
+                #otherwise seed with target
+                else:
+                    seed_model = mdl_start.copy()
+                mdl_lsgst_list = checkpoint.mdl_list
+                final_objfn = checkpoint.final_objfn
+                #final_objfn initialized to None in the GateSetTomographyCheckpoint and will be overwritten
+                #during the loop below unless the last completed iteration is the final iteration
+                #in which case the loop should be skipped. If so I think it is ok that this gets
+                #left set to None. There looks to be some logic for handling this and it looks
+                #like the serialization routines effectively do this already, as the value
+                #of this is lost between writing and reading.
+            else:
+                NotImplementedError('The only currently valid checkpoint inputs are None and GateSetTomographyCheckpoint.')
+            
+            #note the last_completed_iter value is initialized to -1 so the below line 
+            # will have us correctly starting at 0 if this is a fresh checkpoint.
+            starting_idx = checkpoint.last_completed_iter + 1
 
-        tnxt = _time.time(); profiler.add_time('GST: Prep Initial seed', tref); tref = tnxt
+        else:
+            seed_model = mdl_start.copy()
+            mdl_lsgst_list = []
+            starting_idx = 0
 
+        tnxt = _time.time(); profiler.add_time('GST: Prep Initial seed', tref); tref = tnxt
+        
         #Run Long-sequence GST on data
-        mdl_lsgst_list, optimums_list, final_objfn = _alg.run_iterative_gst(
-            ds, mdl_start, bulk_circuit_lists, self.optimizer,
+        #Use the generator based version and query each of the intermediate results.
+        gst_iter_generator = _alg.iterative_gst_generator( 
+            ds, seed_model, bulk_circuit_lists, self.optimizer,
             self.objfn_builders.iteration_builders, self.objfn_builders.final_builders,
-            resource_alloc, printer)
+            resource_alloc, starting_idx, printer)
+
+        #The optima don't actually get used right now, so don't bother trying to
+        #checkpoint these.
+        optima_list = []
+
+        #Now loop through the generator and query the intermediate results:
+        #Do all but the last circuit list.
+        for i in range(starting_idx, len(bulk_circuit_lists)):
+            #then do the final iteration slightly differently since the generator should
+            #give three return values.
+            if i==len(bulk_circuit_lists)-1:
+                mdl_iter, opt_iter, final_objfn =  next(gst_iter_generator)
+            else:
+                mdl_iter, opt_iter =  next(gst_iter_generator)
+            mdl_lsgst_list.append(mdl_iter)
+            optima_list.append(opt_iter)
+
+            if not disable_checkpointing:
+                #update the checkpoint along the way:
+                checkpoint.mdl_list = mdl_lsgst_list
+                checkpoint.last_completed_iter += 1
+                checkpoint.last_completed_circuit_list = bulk_circuit_lists[i]
+                #write the updated checkpoint to disk:
+                if resource_alloc.comm_rank == 0:
+                    checkpoint.write(f'{checkpoint_path}_iteration_{i}.json')
 
         tnxt = _time.time(); profiler.add_time('GST: total iterative optimization', tref); tref = tnxt
-
+    
         #set parameters
         parameters = _collections.OrderedDict()
         parameters['protocol'] = self  # Estimates can hold sub-Protocols <=> sub-results
@@ -1337,7 +1430,7 @@ def run(self, data, memlimit=None, comm=None):
         # Note: we associate 'final_cache' with the Estimate, which means we assume that *all*
         # of the models in the estimate can use same evaltree, have the same default prep/POVMs, etc.
 
-        #TODO: add qtys about fit from optimums_list
+        #TODO: add qtys about fit from optima_list
 
         ret = ModelEstimateResults(data, self)
 
@@ -1581,11 +1674,20 @@ class StandardGST(_proto.Protocol):
         be used.
     """
 
-    def __init__(self, modes="full TP,CPTP,Target", gaugeopt_suite='stdgaugeopt', target_model=None,
+    def __init__(self, modes=('full TP','CPTPLND','Target'), gaugeopt_suite='stdgaugeopt', target_model=None,
                  models_to_test=None, objfn_builders=None, optimizer=None, badfit_options=None, verbosity=2, name=None):
 
         super().__init__(name)
-        self.modes = modes.split(',')
+        if isinstance(modes, str):
+            if ',' in modes:
+                self.modes = modes.split(',')
+                _warnings.warn("The use of a comma-separated string as input for 'modes' is deprecated " 
+                            + " and may be removed in a future release. Please pass in a list or tuple"
+                                +" (or other iterable) of strings")
+            else:
+                self.modes = [modes] #Cast to a list for uniformity
+        else:
+            self.modes = modes
         self.models_to_test = models_to_test
         self.target_model = target_model
         self.gaugeopt_suite = GSTGaugeOptSuite.cast(gaugeopt_suite)
@@ -1614,7 +1716,7 @@ def __init__(self, modes="full TP,CPTP,Target", gaugeopt_suite='stdgaugeopt', ta
     #    data = _proto.ProtocolData(design, dataset)
     #    return self.run(data)
 
-    def run(self, data, memlimit=None, comm=None):
+    def run(self, data, memlimit=None, comm=None, checkpoint= None, checkpoint_path=None, disable_checkpointing = False):
         """
         Run this protocol on `data`.
 
@@ -1630,6 +1732,23 @@ def run(self, data, memlimit=None, comm=None):
             When not ``None``, an MPI communicator used to run this protocol
             in parallel.
 
+        checkpoint : StandardGSTCheckpoint, optional (default None)
+            If specified use a previously generated checkpoint object to restart
+            or warm start this run part way through.
+
+        checkpoint_path : str, optional (default None)
+            A string for the path/name to use for writing intermediate checkpoint
+            files to disk. Format is {path}/{name}, without inclusion of the json
+            file extension. This {path}/{name} combination will have the latest
+            completed iteration number appended to it before writing it to disk.
+            If none, the value of {name} will be set to the name of the protocol
+            being run.
+        
+        disable_checkpointing : bool, optional (default False)
+            When set to True checkpoint objects will not be constructed and written
+            to disk during the course of this protocol. It is strongly recommended
+            that this be kept set to False without good reason to disable the checkpoints.
+
         Returns
         -------
         ProtocolResults
@@ -1659,26 +1778,66 @@ def run(self, data, memlimit=None, comm=None):
         else:
             target_model = None  # Usually this path leads to an error being raised below.
 
+        if not disable_checkpointing:
+            #Set the checkpoint_path variable if None
+            if checkpoint_path is None:
+                checkpoint_path_base = _pathlib.Path('./standard_gst_checkpoints/' + self.name)
+            else:
+                #cast this to a pathlib path with the file extension (suffix) dropped
+                checkpoint_path_base = _pathlib.Path(checkpoint_path).with_suffix('')
+
+            #If there is no checkpoint we should start from with the seed model,
+            #otherwise we should seed the next iteration with the last iteration's result.
+            #If there is no checkpoint initialize mdl_lsgst_list and final_objfn to be empty, 
+            #otherwise re-initialize their values from the checkpoint   
+            if checkpoint is None:
+                checkpoint = StandardGSTCheckpoint(modes)
+                #pre-populate the children of the StandardGSTCheckpoint
+                #with properly initialized checkpoint objects for each of
+                #the protocols.
+                child_dict = {}
+                for mode in modes:
+                    if mode == "Target" or mode in models_to_test:
+                        child_dict[mode] = _ModelTestCheckpoint(name = mode, parent = checkpoint)
+                    else:
+                        child_dict[mode] = GateSetTomographyCheckpoint(name = mode, parent = checkpoint)
+                checkpoint.children = child_dict
+            elif isinstance(checkpoint, StandardGSTCheckpoint):
+                pass
+            else:
+                NotImplementedError('The only currently valid checkpoint inputs are None and StandardGSTCheckpoint.')
+
         ret = ModelEstimateResults(data, self)
         with printer.progress_logging(1):
             for i, mode in enumerate(modes):
                 printer.show_progress(i, len(modes), prefix='-- Std Practice: ', suffix=' (%s) --' % mode)
-
+                if not disable_checkpointing:
+                    #pre python 3.9 compatible version.
+                    checkpoint_path = checkpoint_path_base.with_name(f"{checkpoint_path_base.stem}_{mode.replace(' ', '_')}")
+                    #The line below only works for python 3.9+
+                    #checkpoint_path = checkpoint_path_base.with_stem(f"{checkpoint_path_base.stem}_{mode.replace(' ', '_')}")
                 if mode == "Target":
                     if target_model is None:
                         raise ValueError(("Must specify `target_model` when creating this StandardGST, since one could"
                                           " not be inferred from the given experiment design."))
-
-                    model_to_test = target_model
-                    mdltest = _ModelTest(model_to_test, target_model, self.gaugeopt_suite,
+                    
+                    mdltest = _ModelTest(target_model, target_model, self.gaugeopt_suite,
                                          mt_builder, self.badfit_options, verbosity=printer - 1, name=mode)
-                    result = mdltest.run(data, memlimit, comm)
+                    if not disable_checkpointing:
+                        result = mdltest.run(data, memlimit, comm, checkpoint = checkpoint.children[mode],
+                                            checkpoint_path=checkpoint_path)
+                    else:
+                        result = mdltest.run(data, memlimit, comm, disable_checkpointing=True)
                     ret.add_estimates(result)
 
                 elif mode in models_to_test:
                     mdltest = _ModelTest(models_to_test[mode], target_model, self.gaugeopt_suite,
                                          None, self.badfit_options, verbosity=printer - 1, name=mode)
-                    result = mdltest.run(data, memlimit, comm)
+                    if not disable_checkpointing:
+                        result = mdltest.run(data, memlimit, comm, checkpoint = checkpoint.children[mode],
+                                            checkpoint_path=checkpoint_path)
+                    else:
+                        result = mdltest.run(data, memlimit, comm, disable_checkpointing=True)
                     ret.add_estimates(result)
 
                 else:
@@ -1699,7 +1858,11 @@ def run(self, data, memlimit=None, comm=None):
                     initial_model = GSTInitialModel(initial_model, self.starting_point.get(mode, None))
                     gst = GST(initial_model, self.gaugeopt_suite, self.objfn_builders,
                               self.optimizer, self.badfit_options, verbosity=printer - 1, name=mode)
-                    result = gst.run(data, memlimit, comm)
+                    if not disable_checkpointing:
+                        result = gst.run(data, memlimit, comm, checkpoint = checkpoint.children[mode],
+                                        checkpoint_path=checkpoint_path)
+                    else:
+                        result = gst.run(data, memlimit, comm, disable_checkpointing=True)
                     ret.add_estimates(result)
 
         return ret
@@ -2936,6 +3099,156 @@ def __str__(self):
         s += "  " + "\n  ".join(list(self.estimates.keys())) + "\n"
         s += "\n"
         return s
+    
+class GateSetTomographyCheckpoint(_proto.ProtocolCheckpoint):
+    """
+    A class for storing intermediate results associated with running
+    a GateSetTomography protocol's run method to allow for restarting
+    that method partway through.
+
+    Parameters
+    ----------
+    mdl_list : list of models, optional (default None)
+        Current list of models for each of the completed iterations of the protocol.
+
+    last_completed_iter : int, optional (default -1)
+        Index of the last iteration what was successfully completed.
+
+    last_completed_circuit_list : list of Circuit objects, CircuitList or equivalent, optional (default None)
+        A list of Circuit objects corresponding to the last iteration successfully completed.
+    
+    final_objfn : ModelDatasetCircuitStore, optional (Default None)
+        A ModelDatasetCircuitStore object corresponding to the final evaluated objective function.
+        Not currently serialized or used during the warmstarting, so purely informational and may
+        not always be initialized.
+
+    name : str, optional (default None)
+        An optional name for the checkpoint. Note this is not necessarily the name used in the
+        automatic generation of filenames when written to disk.
+
+    parent : ProtocolCheckpoint, optional (default None)
+        When specified this checkpoint object is treated as the child of another ProtocolCheckpoint
+        object that acts as the parent. When present, the parent's `write` method supersedes
+        the child objects and is called when calling `write` on the child. Currently only used
+        in the implementation of StandardGSTCheckpoint.
+    """
+
+    def __init__(self, mdl_list = None, last_completed_iter = -1, 
+                 last_completed_circuit_list = None, final_objfn = None,
+                 name= None, parent = None):
+
+        self.mdl_list = mdl_list if mdl_list is not None else []
+        self.last_completed_iter = last_completed_iter
+        self.last_completed_circuit_list = last_completed_circuit_list if last_completed_circuit_list is not None else []
+        self.final_objfn = final_objfn
+
+        super().__init__(name, parent)
+
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+        state.update({'mdl_list': [mdl.to_nice_serialization() for mdl in self.mdl_list],
+                      'last_completed_iter': self.last_completed_iter,
+                      'last_completed_circuit_list': [ckt.str for ckt in self.last_completed_circuit_list],
+                      'final_objfn': self.final_objfn,
+                      'name': self.name
+                      })
+        return state
+
+    @classmethod
+    def _from_nice_serialization(cls, state):  # memo holds already de-serialized objects
+        mdl_list = [_Model.from_nice_serialization(mdl) for mdl in state['mdl_list']]
+        last_completed_iter = state['last_completed_iter']
+        last_completed_circuit_list = [Circuit(ckt_str) for ckt_str in state['last_completed_circuit_list']]
+        final_objfn = state['final_objfn']
+        name = state['name']
+        return cls(mdl_list, last_completed_iter, last_completed_circuit_list, 
+                   final_objfn, name)
+    
+
+
+class StandardGSTCheckpoint(_proto.ProtocolCheckpoint):
+    """
+    A class for storing intermediate results associated with running
+    a StandardGST protocol's run method to allow for restarting
+    that method partway through. This class acts as a container
+    class for some set of child GateSetTomographyCheckpoint and
+    ModelTestCheckpoint objects for each of the sub-protocols run in
+    the course of the StandardGST protocol.
+
+    Parameters
+    ----------
+    modes : list of str, optional (default None)
+        A list of strings corresponding to the mode labels being run in the
+        StandardGST protocol object that generated this checkpoint.
+
+    children : dict, optional (default None)
+        A dictionary whose keys correspond to modes (i.e. the same elements as the
+        modes kwarg) and whose values are either GateSetTomographyCheckpoint
+        or ModelTestCheckpoint objects, depending on whichever is appropriate
+        for that mode.
+        
+    name : str, optional (default None)
+        An optional name for the checkpoint. Note this is not necessarily the name used in the
+        automatic generation of filenames when written to disk. 
+    """
+
+    def __init__(self, modes= None, children = None, name= None):
+        self.modes = modes
+        self.children = children if children is not None else {}
+        super().__init__(name)
+    
+    @property
+    def children(self):
+        return self._children
+    
+    @children.setter
+    def children(self, child_dict):
+        self._children = child_dict
+        #also initialize something for storing child types for use in
+        #hacky deserialization
+        self.child_types = {}
+        for mode, child in child_dict.items():
+            if isinstance(child, _ModelTestCheckpoint):
+                self.child_types[mode] = 'modeltest'
+            elif isinstance(child, GateSetTomographyCheckpoint):
+                self.child_types[mode] = 'gatesettomography'
+            else:
+                raise ValueError('StandardGSTCheckpoint currently only supports'\
+                                 +' child checkpoint types that are ModelTestCheckpoint'\
+                                 +' or GateSetTomographyCheckpoint objects.')
+            
+
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+        state.update({'modes': self.modes,
+                      'children': {mode: self.children[mode]._to_nice_serialization() for mode in self.modes},
+                      'child_types': self.child_types,
+                      'name': self.name
+                      })
+        return state
+
+    @classmethod
+    def _from_nice_serialization(cls, state):  # memo holds already de-serialized objects
+        modes = state['modes']
+        child_types = state['child_types']
+        #reinitialize the checkpoint objects for the children
+        child_serializations = state['children']
+        children = {}
+        for mode in modes:
+            if child_types[mode] == 'modeltest':
+                children[mode] = _ModelTestCheckpoint._from_nice_serialization(child_serializations[mode])
+            elif child_types[mode] == 'gatesettomography':
+                children[mode] = GateSetTomographyCheckpoint._from_nice_serialization(child_serializations[mode])
+            else:
+                raise ValueError('StandardGSTCheckpoint currently only supports'\
+                                 +' child checkpoint types that are ModelTestCheckpoint'\
+                                 +' or GateSetTomographyCheckpoint objects.')
+        name = state['name']
+        ret = cls(modes, children, name)
+        #relink the parents for each of the children so they point to the newly returned instance ret
+        for mode in modes:
+            ret.children[mode].parent = ret
+        return ret
 
 
 GSTDesign = GateSetTomographyDesign
diff --git a/pygsti/protocols/protocol.py b/pygsti/protocols/protocol.py
index cf7d7cc63..ae326b23c 100644
--- a/pygsti/protocols/protocol.py
+++ b/pygsti/protocols/protocol.py
@@ -24,6 +24,7 @@
 from pygsti.tools import listtools as _lt
 from pygsti.tools.dataframetools import _process_dataframe
 from pygsti.baseobjs.mongoserializable import MongoSerializable as _MongoSerializable
+from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 
 
 class Protocol(_MongoSerializable):
@@ -3208,6 +3209,36 @@ def run(self, edesign, memlimit=None, comm=None):
                             self.record_zero_counts, comm, memlimit, self.times)
         return ProtocolData(edesign, ds)
 
+class ProtocolCheckpoint(_NicelySerializable):
+    """
+    Class for storing checkpointing intermediate progress during
+    the running of a protocol in order to enable restarting subsequent
+    runs of the protocol from that point.
+    
+    Parameters
+    ----------
+    name : str
+        Name of the protocol associated with this checkpoint.
+
+    parent : ProtocolCheckpoint, optional (default None)
+        When specified this checkpoint object is treated as the child of another ProtocolCheckpoint
+        object that acts as the parent. When present, the parent's `write` method supersedes
+        the child objects and is called when calling `write` on the child. Currently only used
+        in the implementation of StandardGSTCheckpoint.
+    """
+
+    def __init__(self, name, parent = None):
+        self.name = name
+        self.parent = parent
+        #Need to add this for MongoDB serialization related reasons.
+        self._dbcoordinates = None
+
+    def write(self, path):
+        if self.parent is not None:
+            self.parent.write(path)
+        else:
+            super().write(path)
+
 
 #In the future, we could put this function into a base class for
 # the classes that utilize it above, so it would become a proper method.


Test patch:
diff --git a/pygsti/protocols/modeltest.py b/pygsti/protocols/modeltest.py
index 484dbf9cd..e9b691271 100644
--- a/pygsti/protocols/modeltest.py
+++ b/pygsti/protocols/modeltest.py
@@ -12,7 +12,7 @@
 
 import collections as _collections
 import warnings as _warnings
-
+import pathlib as _pathlib
 from pygsti.baseobjs.profiler import DummyProfiler as _DummyProfiler
 from pygsti.objectivefns.objectivefns import ModelDatasetCircuitsStore as _ModelDatasetCircuitStore
 from pygsti.protocols.estimate import Estimate as _Estimate
@@ -20,6 +20,7 @@
 from pygsti import baseobjs as _baseobjs
 from pygsti import models as _models
 from pygsti.objectivefns import objectivefns as _objfns
+from pygsti.circuits import Circuit
 from pygsti.circuits.circuitlist import CircuitList as _CircuitList
 from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
 
@@ -130,7 +131,7 @@ def __init__(self, model_to_test, target_model=None, gaugeopt_suite=None,
     #    design = _StandardGSTDesign(target_model, prep_fiducials, meas_fiducials, germs, maxLengths)
     #    return self.run(_proto.ProtocolData(design, dataset))
 
-    def run(self, data, memlimit=None, comm=None):
+    def run(self, data, memlimit=None, comm=None, checkpoint=None, checkpoint_path=None, disable_checkpointing= False):
         """
         Run this protocol on `data`.
 
@@ -146,6 +147,23 @@ def run(self, data, memlimit=None, comm=None):
             When not ``None``, an MPI communicator used to run this protocol
             in parallel.
 
+        checkpoint : ModelTestCheckpoint, optional (default None)
+            If specified use a previously generated checkpoint object to restart
+            or warm start this run part way through.
+
+        checkpoint_path : str, optional (default None)
+            A string for the path/name to use for writing intermediate checkpoint
+            files to disk. Format is {path}/{name}, without inclusion of the json
+            file extension. This {path}/{name} combination will have the latest
+            completed iteration number appended to it before writing it to disk.
+            If none, the value of {name} will be set to the name of the protocol
+            being run.
+        
+        disable_checkpointing : bool, optional (default False)
+            When set to True checkpoint objects will not be constructed and written
+            to disk during the course of this protocol. It is strongly recommended
+            that this be kept set to False without good reason to disable the checkpoints.
+
         Returns
         -------
         ModelEstimateResults
@@ -176,15 +194,58 @@ def run(self, data, memlimit=None, comm=None):
         bulk_circuit_lists = [_CircuitList(lst, aliases, self.circuit_weights)
                               if not isinstance(lst, _CircuitList) else lst
                               for lst in circuit_lists]
-        objfn_vals = []
-        chi2k_distributed_vals = []
+        
+        if not disable_checkpointing:
+            #Set the checkpoint_path variable if None
+            if checkpoint_path is None:
+                checkpoint_path = _pathlib.Path('./model_test_checkpoints/' + self.name)
+            else:
+                #cast this to a pathlib path with the file extension (suffix) dropped
+                checkpoint_path = _pathlib.Path(checkpoint_path).with_suffix('')
+            
+            #create the parent directory of the checkpoint if needed:
+            checkpoint_path.parent.mkdir(parents=True, exist_ok=True)
+            
+            #If there is no checkpoint we should start from with the first circuit list,
+            #otherwise we should load in the cached results and start from the point
+            # in the objective function calculation we left off.
+            if checkpoint is None:
+                objfn_vals = []
+                chi2k_distributed_vals = []
+                checkpoint = ModelTestCheckpoint()
+            elif isinstance(checkpoint, ModelTestCheckpoint):
+                objfn_vals = checkpoint.objfn_vals
+                chi2k_distributed_vals = checkpoint.chi2k_distributed_vals
+            else:
+                NotImplementedError('The only currently valid checkpoint inputs are None and ModelTestCheckpoint.')
+
+            #Check the last completed iteration identified in the checkpoint and set that as the
+            #starting point for the iteration through bulk_circuit_lists. This starts at -1 for
+            #a freshly initialized checkpoint.
+            starting_idx = checkpoint.last_completed_iter + 1
+        else:
+            starting_idx = 0
+            objfn_vals = []
+            chi2k_distributed_vals = []
+
         assert(len(self.objfn_builders) == 1), "Only support for a single objective function so far."
-        for circuit_list in bulk_circuit_lists:
+        for i in range(starting_idx, len(bulk_circuit_lists)):
+            circuit_list = bulk_circuit_lists[i]
             objective = self.objfn_builders[0].build(the_model, ds, circuit_list, resource_alloc, printer - 1)
             f = objective.fn(the_model.to_vector())
             objfn_vals.append(f)
             chi2k_distributed_vals.append(objective.chi2k_distributed_qty(f))
 
+            if not disable_checkpointing:
+                #Update the checkpoint:
+                checkpoint.objfn_vals = objfn_vals
+                checkpoint.chi2k_distributed_vals = chi2k_distributed_vals
+                checkpoint.last_completed_iter += 1
+                checkpoint.last_completed_circuit_list= circuit_list
+                #write the updated checkpoint to disk:
+                if resource_alloc.comm_rank == 0:
+                    checkpoint.write(f'{checkpoint_path}_iteration_{i}.json')
+
         mdc_store = _ModelDatasetCircuitStore(the_model, ds, bulk_circuit_lists[-1], resource_alloc)
         parameters = _collections.OrderedDict()
         parameters['final_objfn_builder'] = self.objfn_builders[-1]
@@ -211,3 +272,69 @@ def run(self, data, memlimit=None, comm=None):
         return _add_gaugeopt_and_badfit(ret, self.name, target_model, self.gaugeopt_suite,
                                         self.unreliable_ops, self.badfit_options,
                                         None, resource_alloc, printer)
+
+
+class ModelTestCheckpoint(_proto.ProtocolCheckpoint):
+    """
+    A class for storing intermediate results associated with running
+    a ModelTest protocol's run method to allow for restarting
+    that method partway through.
+
+    Parameters
+    ----------
+    last_completed_iter : int, optional (default -1)
+        Index of the last iteration what was successfully completed.
+
+    last_completed_circuit_list : list of Circuit objects, CircuitList or equivalent, optional (default None)
+        A list of Circuit objects corresponding to the last iteration successfully completed.
+
+    objfn_vals : list, optional (default None)
+        A list of the current objective function values for each iteration/circuit list
+        evaluated during the ModelTest protocol.
+
+    chi2k_distributed_vals : list, optional (default None)
+        A list of the current objective function values for each iteration/circuit list
+         evaluated during the ModelTest protocol rescaled so as to have an expected chi-squared
+         distribution under the null hypothesis of Wilks' theorem.
+        
+    name : str, optional (default None)
+        An optional name for the checkpoint. Note this is not necessarily the name used in the
+        automatic generation of filenames when written to disk. 
+    
+    parent : ProtocolCheckpoint, optional (default None)
+        When specified this checkpoint object is treated as the child of another ProtocolCheckpoint
+        object that acts as the parent. When present, the parent's `write` method supersedes
+        the child objects and is called when calling `write` on the child. Currently only used
+        in the implementation of StandardGSTCheckpoint.
+
+    """
+
+    def __init__(self, last_completed_iter = -1, 
+                 last_completed_circuit_list = None, objfn_vals = None,
+                 chi2k_distributed_vals=None, name= None, parent = None):
+        self.last_completed_iter = last_completed_iter
+        self.last_completed_circuit_list = last_completed_circuit_list if last_completed_circuit_list is not None else []
+        self.objfn_vals = objfn_vals if objfn_vals is not None else []
+        self.chi2k_distributed_vals = chi2k_distributed_vals if chi2k_distributed_vals is not None else []
+
+        super().__init__(name, parent)
+
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+        state.update({'last_completed_iter': self.last_completed_iter,
+                      'last_completed_circuit_list': [ckt.str for ckt in self.last_completed_circuit_list],
+                      'objfn_vals': self.objfn_vals,
+                      'chi2k_distributed_vals': self.chi2k_distributed_vals,
+                      'name': self.name
+                      })
+        return state
+
+    @classmethod
+    def _from_nice_serialization(cls, state):  # memo holds already de-serialized objects
+        last_completed_iter = state['last_completed_iter']
+        last_completed_circuit_list = [Circuit(ckt_str) for ckt_str in state['last_completed_circuit_list']]
+        objfn_vals = state['objfn_vals']
+        chi2k_distributed_vals = state['chi2k_distributed_vals']
+        name = state['name']
+        return cls(last_completed_iter, last_completed_circuit_list, 
+                   objfn_vals, chi2k_distributed_vals, name)
diff --git a/test/test_packages/drivers/test_drivers.py b/test/test_packages/drivers/test_drivers.py
index 68a58a309..9b3f8533b 100644
--- a/test/test_packages/drivers/test_drivers.py
+++ b/test/test_packages/drivers/test_drivers.py
@@ -1,12 +1,11 @@
-import unittest
-
+import unittest # noqa: E999
+import numpy as np
 from pygsti.forwardsims.mapforwardsim import MapForwardSimulator
-
+from numpy.linalg import norm
 import pygsti
 from pygsti.modelpacks.legacy import std1Q_XYI as std
 from ..testutils import BaseTestCase, compare_files, temp_files, regenerate_references
 
-
 class DriversTestCase(BaseTestCase):
 
     def setUp(self):
@@ -62,7 +61,8 @@ def test_longSequenceGST_fiducialPairReduction(self):
             std.target_model(), std.fiducials, std.fiducials, std.germs, maxLens,
             fid_pairs=fidPairs)
 
-        result = pygsti.run_long_sequence_gst_base(ds, std.target_model(), gfprStructs, verbosity=0)
+        result = pygsti.run_long_sequence_gst_base(ds, std.target_model(), gfprStructs, verbosity=0,
+                                                   disable_checkpointing = True)
         pygsti.report.create_standard_report(result, temp_files + "/full_report_GFPR",
                                              "GFPR report", verbosity=2)
 
@@ -88,7 +88,8 @@ def test_longSequenceGST_fiducialPairReduction(self):
             std.target_model(), std.fiducials, std.fiducials, std.germs, maxLens,
             fid_pairs=fidPairsDict)
 
-        result = pygsti.run_long_sequence_gst_base(ds, std.target_model(), pfprStructs, verbosity=0)
+        result = pygsti.run_long_sequence_gst_base(ds, std.target_model(), pfprStructs, verbosity=0,
+                                                   disable_checkpointing = True)
         pygsti.report.create_standard_report(result, temp_files + "/full_report_PFPR",
                                              "PFPR report", verbosity=2)
 
@@ -106,7 +107,8 @@ def test_longSequenceGST_randomReduction(self):
             maxLens, fidPairs, ts, keep_fraction=0.5, keep_seed=1234)
         result = self.runSilent(pygsti.run_long_sequence_gst_base,
                                 ds, std.target_model(), reducedLists,
-                                advanced_options={'truncScheme': ts})
+                                advanced_options={'truncScheme': ts},
+                                disable_checkpointing=True)
 
         #create a report...
         pygsti.report.create_standard_report(result, temp_files + "/full_report_RFPR",
@@ -120,7 +122,8 @@ def test_longSequenceGST_randomReduction(self):
             maxLens, fidPairs, ts, keep_fraction=0.5, keep_seed=1234)
         result2 = self.runSilent(pygsti.run_long_sequence_gst_base,
                                  ds, std.target_model(), reducedLists,
-                                 advanced_options={'truncScheme': ts})
+                                 advanced_options={'truncScheme': ts},
+                                 disable_checkpointing=True)
 
         #create a report...
         pygsti.report.create_standard_report(result2, temp_files + "/full_report_RFPR2.html",
@@ -130,12 +133,12 @@ def test_longSequenceGST_CPTP(self):
         ds = pygsti.data.DataSet(file_to_load_from=compare_files + "/drivers.dataset")
 
         target_model = std.target_model()
-        target_model.set_all_parameterizations("CPTP")
+        target_model.set_all_parameterizations("CPTPLND")
 
         maxLens = self.maxLens
         result = self.runSilent(pygsti.run_long_sequence_gst,
                                 ds, target_model, std.fiducials, std.fiducials,
-                                std.germs, maxLens)
+                                std.germs, maxLens, disable_checkpointing=True)
 
         #create a report...
         pygsti.report.create_standard_report(result, temp_files + "/full_report_CPTPGates",
@@ -151,7 +154,7 @@ def test_longSequenceGST_Sonly(self):
         maxLens = self.maxLens
         result = self.runSilent(pygsti.run_long_sequence_gst,
                                 ds, target_model, std.fiducials, std.fiducials,
-                                std.germs, maxLens)
+                                std.germs, maxLens, disable_checkpointing=True)
 
         #create a report...
         pygsti.report.create_standard_report(result, temp_files + "/full_report_SGates.html",
@@ -174,7 +177,7 @@ def test_longSequenceGST_GLND(self):
         maxLens = self.maxLens
         result = self.runSilent(pygsti.run_long_sequence_gst,
                                 ds, target_model, std.fiducials, std.fiducials,
-                                std.germs, maxLens)
+                                std.germs, maxLens, disable_checkpointing=True)
 
         #create a report...
         pygsti.report.create_standard_report(result, temp_files + "/full_report_SGates",
@@ -190,7 +193,7 @@ def test_longSequenceGST_HplusS(self):
         maxLens = self.maxLens
         result = self.runSilent(pygsti.run_long_sequence_gst,
                                 ds, target_model, std.fiducials, std.fiducials,
-                                std.germs, maxLens)
+                                std.germs, maxLens, disable_checkpointing=True)
 
         #create a report...
         pygsti.report.create_standard_report(result, temp_files + "/full_report_HplusSGates",
@@ -205,7 +208,8 @@ def test_longSequenceGST_badfit(self):
         maxLens = self.maxLens
         result = self.runSilent(pygsti.run_long_sequence_gst,
                                 ds, std.target_model(), std.fiducials, std.fiducials,
-                                std.germs, maxLens, advanced_options={'bad_fit_threshold': -100})
+                                std.germs, maxLens, advanced_options={'bad_fit_threshold': -100},
+                                disable_checkpointing=True)
 
         pygsti.report.create_standard_report(result, temp_files + "/full_report_badfit",
                                              "badfit report", verbosity=2)
@@ -218,9 +222,10 @@ def test_stdpracticeGST(self):
         maxLens = self.maxLens
         result = self.runSilent(pygsti.run_stdpractice_gst,
                                 ds, std.target_model().create_processor_spec(), std.fiducials, std.fiducials,
-                                std.germs, maxLens, modes="full TP,CPTP,Test,Target",
+                                std.germs, maxLens, modes=['full TP','CPTPLND','Test','Target'],
                                 models_to_test = {"Test": mdl_guess},
-                                comm=None, mem_limit=None, verbosity=5)
+                                comm=None, mem_limit=None, verbosity=5,
+                                disable_checkpointing=True)
         pygsti.report.create_standard_report(result, temp_files + "/full_report_stdpractice",
                                              "Std Practice Test Report", verbosity=2)
 
@@ -243,6 +248,98 @@ def test_bootstrap(self):
                 2, ds_defaultMaxLens, 'parametric', std.fiducials, std.fiducials,
                 std.germs, default_maxLens, input_model=mdl, target_model=tp_target,
                 return_data=False) #test when max_lengths == None ?
+                
+    def test_GST_checkpointing(self):
+        ds = pygsti.data.DataSet(file_to_load_from=compare_files + "/drivers.dataset")
+        maxLens = self.maxLens
+
+        #Make list-of-lists of GST operation sequences
+        fullStructs = pygsti.circuits.make_lsgst_structs(
+            std.target_model(), std.fiducials, std.fiducials, std.germs, maxLens)
+        
+        #Test GateSetTomographyCheckpoint:
+        #First run from scratch:
+        result_gst = pygsti.run_long_sequence_gst_base(ds, std.target_model(), fullStructs, verbosity=0, 
+                                                   checkpoint_path= temp_files + '/checkpoint_testing/GateSetTomography')
+                                                   
+        #double check that we can read in this checkpoint object correctly:
+        gst_checkpoint = pygsti.protocols.GateSetTomographyCheckpoint.read(temp_files + '/checkpoint_testing/GateSetTomography_iteration_0.json')
+        
+        #run GST using this checkpoint
+        result_gst_warmstart = pygsti.run_long_sequence_gst_base(ds, std.target_model(), fullStructs, verbosity=0,
+                                                                 checkpoint = gst_checkpoint,
+                                                                 checkpoint_path= temp_files + '/checkpoint_testing/GateSetTomography')
+        
+        diff = norm(result_gst.estimates['GateSetTomography'].models['final iteration estimate'].to_vector()-
+                         result_gst_warmstart.estimates['GateSetTomography'].models['final iteration estimate'].to_vector())
+        #Assert that this gives the same result as before:
+        self.assertTrue(diff<=1e-10)
+            
+
+    def test_ModelTest_checkpointing(self):
+        ds = pygsti.data.DataSet(file_to_load_from=compare_files + "/drivers.dataset")
+        maxLens = self.maxLens
+
+        #Next test ModelTestCheckpoint
+        #First run from scratch:
+        result_modeltest = pygsti.run_model_test(std.target_model(), ds,std.target_model().create_processor_spec(), 
+                                                 std.fiducials, std.fiducials, std.germs,
+                                                 maxLens, verbosity=0, 
+                                                 checkpoint_path= temp_files + '/checkpoint_testing/ModelTest')
+                                                   
+        #double check that we can read in this checkpoint object correctly:
+        model_test_checkpoint = pygsti.protocols.ModelTestCheckpoint.read(temp_files + '/checkpoint_testing/ModelTest_iteration_0.json')
+        
+        #run GST using this checkpoint
+        result_modeltest_warmstart = pygsti.run_model_test(std.target_model(), ds,std.target_model().create_processor_spec(), 
+                                                           std.fiducials, std.fiducials, std.germs,
+                                                           maxLens, verbosity=0,
+                                                           checkpoint = model_test_checkpoint,
+                                                           checkpoint_path= temp_files + '/checkpoint_testing/ModelTest')
+        
+        diff = norm(np.array(result_modeltest.estimates['ModelTest'].parameters['model_test_values'])- 
+                         np.array(result_modeltest_warmstart.estimates['ModelTest'].parameters['model_test_values']))
+        #Assert that this gives the same result as before:
+        self.assertTrue(diff<=1e-10)
+                      
+
+    
+    def test_StandardGST_checkpointing(self):
+        ds = pygsti.data.DataSet(file_to_load_from=compare_files + "/drivers.dataset")
+        maxLens = self.maxLens
+
+        #Finally test StandardGSTCheckpoint
+        #First run from scratch:
+        mdl_guess = std.target_model().depolarize(op_noise=0.01,spam_noise=0.01)
+                
+        result_standardgst = pygsti.run_stdpractice_gst(ds, std.target_model().create_processor_spec(), std.fiducials, std.fiducials,
+                                                        std.germs, maxLens, modes=['full TP','CPTPLND','Test','Target'],
+                                                        models_to_test = {"Test": mdl_guess},
+                                                        comm=None, mem_limit=None, verbosity=0,
+                                                        checkpoint_path= temp_files + '/checkpoint_testing/StandardGST')
+                                                   
+        #double check that we can read in this checkpoint object correctly:
+        standardgst_checkpoint = pygsti.protocols.StandardGSTCheckpoint.read(temp_files + '/checkpoint_testing/StandardGST_CPTPLND_iteration_1.json')
+        
+        #run GST using this checkpoint
+        result_standardgst_warmstart = pygsti.run_stdpractice_gst(ds, std.target_model().create_processor_spec(), std.fiducials, std.fiducials,
+                                                                  std.germs, maxLens, modes=['full TP','CPTPLND','Test','Target'],
+                                                                  models_to_test = {"Test": mdl_guess},
+                                                                  comm=None, mem_limit=None, verbosity=0,
+                                                                  checkpoint = standardgst_checkpoint,
+                                                                  checkpoint_path= temp_files + '/checkpoint_testing/StandardGST')
+
+        #Assert that this gives the same result as before:
+        #diff = norm(result_standardgst.estimates['CPTPLND'].models['final iteration estimate'].to_vector()- 
+        #                 result_standardgst_warmstart.estimates['CPTPLND'].models['final iteration estimate'].to_vector())
+        diff = pygsti.tools.logl(result_standardgst.estimates['CPTPLND'].models['final iteration estimate'], ds)- \
+               pygsti.tools.logl(result_standardgst_warmstart.estimates['CPTPLND'].models['final iteration estimate'], ds)
+               
+        diff1 = norm(result_standardgst.estimates['full TP'].models['final iteration estimate'].to_vector()- 
+                     result_standardgst_warmstart.estimates['full TP'].models['final iteration estimate'].to_vector())
+        
+        self.assertTrue(abs(diff)<=1e-6)
+        self.assertTrue(diff1<=1e-10)
 
 
 if __name__ == "__main__":
diff --git a/test/test_packages/objects/test_hessian.py b/test/test_packages/objects/test_hessian.py
index 1f6750dc9..0bdb48f7f 100644
--- a/test/test_packages/objects/test_hessian.py
+++ b/test/test_packages/objects/test_hessian.py
@@ -131,7 +131,7 @@ def test_confidenceRegion(self):
         edesign = proto.CircuitListsDesign([pygsti.circuits.CircuitList(circuit_struct)
                                             for circuit_struct in self.gss])
         data = proto.ProtocolData(edesign, self.ds)
-        res = proto.ModelEstimateResults(data, proto.StandardGST(modes="TP"))
+        res = proto.ModelEstimateResults(data, proto.StandardGST(modes="full TP"))
 
         #Add estimate for hessian-based CI --------------------------------------------------
         builder = pygsti.objectivefns.PoissonPicDeltaLogLFunction.builder()
@@ -366,7 +366,7 @@ def test_pickle_ConfidenceRegion(self):
         edesign = proto.CircuitListsDesign([pygsti.circuits.CircuitList(circuit_struct)
                                             for circuit_struct in self.gss])
         data = proto.ProtocolData(edesign, self.ds)
-        res = proto.ModelEstimateResults(data, proto.StandardGST(modes="TP"))
+        res = proto.ModelEstimateResults(data, proto.StandardGST(modes="full TP"))
 
         res.add_estimate(
             proto.estimate.Estimate.create_gst_estimate(
diff --git a/test/unit/algorithms/test_core.py b/test/unit/algorithms/test_core.py
index 456617aab..8ccb41a83 100644
--- a/test/unit/algorithms/test_core.py
+++ b/test/unit/algorithms/test_core.py
@@ -410,3 +410,38 @@ def test_do_mlgst_base_forcefn_grad(self):
                                         resource_alloc=None
                                         )
         # TODO assert correctness
+
+    #Add a test for the new generator method.
+    #run_iterative_gst uses this under the hood, so we really only need to separately
+    #test starting at a different starting index in the fits.
+
+    def test_iterative_gst_generator_starting_index(self):
+        generator = core.iterative_gst_generator(
+            self.ds, self.mdl_clgst, self.lsgstStrings,
+            optimizer={'tol': 1e-5},
+            iteration_objfn_builders=['chi2'],
+            final_objfn_builders=['logl'],
+            resource_alloc=None, starting_index=0, verbosity=0
+        )
+
+        models = []
+
+        for _ in range(len(self.lsgstStrings)):
+            models.append(next(generator)[0])
+
+        #now make a new generator starting from a different index
+        #with a start model corresponding to the model in models for
+        #that index.
+        generator1 = core.iterative_gst_generator(
+            self.ds, models[0], self.lsgstStrings,
+            optimizer={'tol': 1e-5},
+            iteration_objfn_builders=['chi2'],
+            final_objfn_builders=['logl'],
+            resource_alloc=None, starting_index=1, verbosity=0
+        )
+        models1= []
+        for _ in range(1,len(self.lsgstStrings)):
+            models1.append(next(generator1)[0])
+
+        #Make sure we get the same result in both cases.
+        self.assertArraysAlmostEqual(models[-1].to_vector(), models1[-1].to_vector())
diff --git a/test/unit/protocols/test_gst.py b/test/unit/protocols/test_gst.py
index 2e56552c6..9ac5ddea2 100644
--- a/test/unit/protocols/test_gst.py
+++ b/test/unit/protocols/test_gst.py
@@ -256,7 +256,7 @@ class StandardGSTTester(BaseProtocolData, BaseCase):
     """
 
     def test_run(self):
-        proto = gst.StandardGST(modes="full TP,CPTPLND,Target")
+        proto = gst.StandardGST(modes=["full TP","CPTPLND","Target"])
         results = proto.run(self.gst_data)
 
         mdl_result = results.estimates["full TP"].models['stdgaugeopt']