Repository: sandialabs/pyGSTi (language: python)

Issue:
Faster Circuit Primitives
I have always been slightly perturbed by the amount of time that things such as COPA layout creation take when performing GST analysis. There are a number of reasons for this, but one family of reasons is that layout generation relies heavily on a number of Circuit primitives that are not as performant as they can be. This PR is a step toward rectifying this through re-implementation of a number of methods identified as bottlenecks through profiling in the context of actual end-to-end GST analysis. 

In terms of the top-level numbers, for 1-qubit GST fitting using the 'full TP' parameterization with an experiment design containing 1120 circuits (max depth 128), the changes on this branch collectively reduce the time required for layout creation by ~50%, and reduce the overall runtime of the end-to-end protocol by ~15%. Your mileage will of course vary, however, and I would anticipate higher relative savings for larger/deeper experiment designs. I have not done profiling for the two-qubit case yet, so I won't posit any guesses on the impact there (I would guess it to be less relatively speaking, but can't easily ballpark this without actually testing).

### Summary of Changes:
Here is a high-level summary of the changes. I will also be leaving some comments on particular bits of code to explain any more detailed reasoning for some changes (when in doubt the answer is likely that I tried a few different things and picked the fastest way to do something).

**Note on synthetic benchmarks:** For the quoted figures below these were mostly done using a randomly generated set of test circuits generated according to the following code snippet:

```
import pygsti
from pygsti.circuits.circuitconstruction import list_random_circuits_onelen
from pygsti.modelpacks import smq1Q_XYI
import numpy as np
from numpy.random import default_rng
from itertools import product
lens = [1,2,4,8,16,32,64]
op_labels = tuple(smq1Q_XYI.target_model().operations.keys())
random_circuits = []
for length in lens:
    random_circuits.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 1234))
random_circuits_2 = []
for length in lens:
    random_circuits_2.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 12345))
rng = default_rng(1234)
rng.shuffle(random_circuits)
rng = default_rng(12345)
rng.shuffle(random_circuits_2)
random_pairs = list(product(random_circuits,repeat=2))
```

For equality testing this was run against `random_pairs`, for inclusion testing against `random_circuits` and `random_circuits_2`, and for other tests against `random_circuits` unless stated otherwise.

1. Circuit equality testing: The primary change was adding a length check to allow quicker short circuiting. A somewhat faster codepath was also added for comparisons between static circuits, avoiding recomputation of a circuit's tuple representation. In synthetic testing this gave a ~3X speed up.
2. Inclusion testing:  Nothing specific was added here, but the change to circuit equality testing also sped this up by a slightly larger factor of ~3.4.
3. Hashing: Hashes are now precomputed and cached at the point at which a circuit is made static (this is an accepted practice in python for immutable objects, e.g. strings cache their hashes). Related to this change, at staticization time we precompute the tuple representation of the circuit and store this as well. This is the thing which actually gets hashed. Previously the `tup` property was called whenever this was needed and this property call regenenerated the tuple representation from scratch each time, which is expensive. In synthetic tests this resulted in a ~4.5X speedup, but in the context of the GST profiling described above this was closer to 10X (the computational savings should scale with circuit depth, so this isn't too shocking of a difference). As part of this change I have added static assertions to a few methods previously missing them for safety (we don't want to inadvertently invalidate the cached hash).
4. Copying: Previously copying a circuit went through the standard `__init__` pathway to create a copied circuit instance. With the new implementation we utilize a much faster specialized codepath inspired by the previously existing `_fastinit` one. But with a few more microptimizations to take advantage of some of the additional promises we can rely on in the context of copying (which we can't as easily in `_fastinit`). Speedups here range from 3-15X depending on the input and output circuit types. With the low end corresponding to switching static to editable (or vice versa), and the high end keeping the input type unchanged.

This is the end of the changes where I profiled these using synthetic tests. For the rest I did this in the context of GST optimization using cProfile/snakeviz, but didn't carefully record before/after relative changes.

5. Microoptimization galore: There are a bunch of microoptimizations of the sort that usually don't matter except for when  code is particularly hot, and the circuit code is. I won't give an exhaustive list, but this includes things like:
    - Attribute access is faster than method/property access/resolution (by 2X or more in some cases). So many trivial (and some non-trivial) property calls have had the corresponding attributes inlined.
    - Inlining functions to reduce function call overhead.
    - List comprehensions are faster than generator expressions. So, e.g. `tuple([i for i in iterable])` is measurably faster than `tuple(i for i in iterable)` in practice.
    - List comprehensions are faster than the 'build via appends in for loop' approach. Often by a notable amount.
    - Using falsy/truthiness directly is faster than the `if len(obj)>0` idiom if all you're doing is checking if an iterable is empty.

6. Rework SeparatePOVMCircuit internals: The `full_effect_labels` property was getting called a lot, and each time it was it was regenerating the result from scratch (with some non-trivial preprocessing). I've updated the internals to cache this, and also deal with updating the cached value when necessary. This gave a measurable speedup.
7. New `sandwich` method for circuits. This takes as input two tuples of Labels which are prepended and appended to a circuit. This is around twice as fast as doing this with a pair of calls to `__add__` (since we only need 1 call to `_fastinit`). This is not an uncommon task (I added as part of some changes to the circuit completion code in OpModel), so if you spot any places where this is being done it'd be worth switching (maybe @rileyjmurray can write a fancy regex to find these instances).

Additional changes (either outside of the Circuit class, or otherwise not directly performance focused):

8. Bulk versions of `complete_circuit` and `split_circuit`: There are now versions of the `complete_circuit` and `split_circuit` methods for the OpModel class, which respectively add prep and povm labels to a circuit or split them off, which operate on a complete list of circuits all at once. This allows for avoiding some model specific lookups (for values which don't change circuit-to-circuit) many times and is faster overall.
9. `expand_subcircuits_inplace`: The implementation of `expand_subcircuits_inplace` has been reworked so that the behavior is now to recursively expand CircuitLabel objects. Previously the expansion routine only took a single pass through the circuit, meaning that if the circuit contained any nested CircuitLabels (which are perfectly happy with being nested, an easy way to get something like this is to call `to_label` on a circuit containing CircuitLabels) we would only end up expanding out the top level. I noticed there were a few places in the layout creation/forward simulation code that implicitly assumed only a single-level of CircuitLabels would be present, which would likely have eventually led to some really annoying bugs to track down. This change covers that edge case, and should only be slightly more expensive than the previous routine in the typical case where there is only a single-level of CircuitLabels.
10. I've added a deprecation warning to the `simulate` method. If there is strong opposition to removing this functionality we can discuss this, of course. My reasoning here is that given the rest of pygsti's structure this functionality really isn't something we should be asking a Circuit to do, as the preferred pattern nowadays is delegate this job to the Model/ForwardSimulator, rather than asking the circuit itself. I don't know how widespread the use of this is (my guess is that this is largely vestigial), but I'd recommend slating this for removal in 0.9.14. Related note, I don't know if we have a running list of things we've deprecated and their expiry date, but if not we should assemble one.

Finally, there are a number of miscellaneous documentation fixes, bugfixes, and updates to unit tests.

P.S. There are some more layout creation performance enhancements inbound on a separate branch, but I've opted to keep these changes together in a separate PR to keep things thematic and make it easier to review without context switching.
Faster Circuit Primitives
I have always been slightly perturbed by the amount of time that things such as COPA layout creation take when performing GST analysis. There are a number of reasons for this, but one family of reasons is that layout generation relies heavily on a number of Circuit primitives that are not as performant as they can be. This PR is a step toward rectifying this through re-implementation of a number of methods identified as bottlenecks through profiling in the context of actual end-to-end GST analysis. 

In terms of the top-level numbers, for 1-qubit GST fitting using the 'full TP' parameterization with an experiment design containing 1120 circuits (max depth 128), the changes on this branch collectively reduce the time required for layout creation by ~50%, and reduce the overall runtime of the end-to-end protocol by ~15%. Your mileage will of course vary, however, and I would anticipate higher relative savings for larger/deeper experiment designs. I have not done profiling for the two-qubit case yet, so I won't posit any guesses on the impact there (I would guess it to be less relatively speaking, but can't easily ballpark this without actually testing).

### Summary of Changes:
Here is a high-level summary of the changes. I will also be leaving some comments on particular bits of code to explain any more detailed reasoning for some changes (when in doubt the answer is likely that I tried a few different things and picked the fastest way to do something).

**Note on synthetic benchmarks:** For the quoted figures below these were mostly done using a randomly generated set of test circuits generated according to the following code snippet:

```
import pygsti
from pygsti.circuits.circuitconstruction import list_random_circuits_onelen
from pygsti.modelpacks import smq1Q_XYI
import numpy as np
from numpy.random import default_rng
from itertools import product
lens = [1,2,4,8,16,32,64]
op_labels = tuple(smq1Q_XYI.target_model().operations.keys())
random_circuits = []
for length in lens:
    random_circuits.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 1234))
random_circuits_2 = []
for length in lens:
    random_circuits_2.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 12345))
rng = default_rng(1234)
rng.shuffle(random_circuits)
rng = default_rng(12345)
rng.shuffle(random_circuits_2)
random_pairs = list(product(random_circuits,repeat=2))
```

For equality testing this was run against `random_pairs`, for inclusion testing against `random_circuits` and `random_circuits_2`, and for other tests against `random_circuits` unless stated otherwise.

1. Circuit equality testing: The primary change was adding a length check to allow quicker short circuiting. A somewhat faster codepath was also added for comparisons between static circuits, avoiding recomputation of a circuit's tuple representation. In synthetic testing this gave a ~3X speed up.
2. Inclusion testing:  Nothing specific was added here, but the change to circuit equality testing also sped this up by a slightly larger factor of ~3.4.
3. Hashing: Hashes are now precomputed and cached at the point at which a circuit is made static (this is an accepted practice in python for immutable objects, e.g. strings cache their hashes). Related to this change, at staticization time we precompute the tuple representation of the circuit and store this as well. This is the thing which actually gets hashed. Previously the `tup` property was called whenever this was needed and this property call regenenerated the tuple representation from scratch each time, which is expensive. In synthetic tests this resulted in a ~4.5X speedup, but in the context of the GST profiling described above this was closer to 10X (the computational savings should scale with circuit depth, so this isn't too shocking of a difference). As part of this change I have added static assertions to a few methods previously missing them for safety (we don't want to inadvertently invalidate the cached hash).
4. Copying: Previously copying a circuit went through the standard `__init__` pathway to create a copied circuit instance. With the new implementation we utilize a much faster specialized codepath inspired by the previously existing `_fastinit` one. But with a few more microptimizations to take advantage of some of the additional promises we can rely on in the context of copying (which we can't as easily in `_fastinit`). Speedups here range from 3-15X depending on the input and output circuit types. With the low end corresponding to switching static to editable (or vice versa), and the high end keeping the input type unchanged.

This is the end of the changes where I profiled these using synthetic tests. For the rest I did this in the context of GST optimization using cProfile/snakeviz, but didn't carefully record before/after relative changes.

5. Microoptimization galore: There are a bunch of microoptimizations of the sort that usually don't matter except for when  code is particularly hot, and the circuit code is. I won't give an exhaustive list, but this includes things like:
    - Attribute access is faster than method/property access/resolution (by 2X or more in some cases). So many trivial (and some non-trivial) property calls have had the corresponding attributes inlined.
    - Inlining functions to reduce function call overhead.
    - List comprehensions are faster than generator expressions. So, e.g. `tuple([i for i in iterable])` is measurably faster than `tuple(i for i in iterable)` in practice.
    - List comprehensions are faster than the 'build via appends in for loop' approach. Often by a notable amount.
    - Using falsy/truthiness directly is faster than the `if len(obj)>0` idiom if all you're doing is checking if an iterable is empty.

6. Rework SeparatePOVMCircuit internals: The `full_effect_labels` property was getting called a lot, and each time it was it was regenerating the result from scratch (with some non-trivial preprocessing). I've updated the internals to cache this, and also deal with updating the cached value when necessary. This gave a measurable speedup.
7. New `sandwich` method for circuits. This takes as input two tuples of Labels which are prepended and appended to a circuit. This is around twice as fast as doing this with a pair of calls to `__add__` (since we only need 1 call to `_fastinit`). This is not an uncommon task (I added as part of some changes to the circuit completion code in OpModel), so if you spot any places where this is being done it'd be worth switching (maybe @rileyjmurray can write a fancy regex to find these instances).

Additional changes (either outside of the Circuit class, or otherwise not directly performance focused):

8. Bulk versions of `complete_circuit` and `split_circuit`: There are now versions of the `complete_circuit` and `split_circuit` methods for the OpModel class, which respectively add prep and povm labels to a circuit or split them off, which operate on a complete list of circuits all at once. This allows for avoiding some model specific lookups (for values which don't change circuit-to-circuit) many times and is faster overall.
9. `expand_subcircuits_inplace`: The implementation of `expand_subcircuits_inplace` has been reworked so that the behavior is now to recursively expand CircuitLabel objects. Previously the expansion routine only took a single pass through the circuit, meaning that if the circuit contained any nested CircuitLabels (which are perfectly happy with being nested, an easy way to get something like this is to call `to_label` on a circuit containing CircuitLabels) we would only end up expanding out the top level. I noticed there were a few places in the layout creation/forward simulation code that implicitly assumed only a single-level of CircuitLabels would be present, which would likely have eventually led to some really annoying bugs to track down. This change covers that edge case, and should only be slightly more expensive than the previous routine in the typical case where there is only a single-level of CircuitLabels.
10. I've added a deprecation warning to the `simulate` method. If there is strong opposition to removing this functionality we can discuss this, of course. My reasoning here is that given the rest of pygsti's structure this functionality really isn't something we should be asking a Circuit to do, as the preferred pattern nowadays is delegate this job to the Model/ForwardSimulator, rather than asking the circuit itself. I don't know how widespread the use of this is (my guess is that this is largely vestigial), but I'd recommend slating this for removal in 0.9.14. Related note, I don't know if we have a running list of things we've deprecated and their expiry date, but if not we should assemble one.

Finally, there are a number of miscellaneous documentation fixes, bugfixes, and updates to unit tests.

P.S. There are some more layout creation performance enhancements inbound on a separate branch, but I've opted to keep these changes together in a separate PR to keep things thematic and make it easier to review without context switching.
Faster Circuit Primitives
I have always been slightly perturbed by the amount of time that things such as COPA layout creation take when performing GST analysis. There are a number of reasons for this, but one family of reasons is that layout generation relies heavily on a number of Circuit primitives that are not as performant as they can be. This PR is a step toward rectifying this through re-implementation of a number of methods identified as bottlenecks through profiling in the context of actual end-to-end GST analysis. 

In terms of the top-level numbers, for 1-qubit GST fitting using the 'full TP' parameterization with an experiment design containing 1120 circuits (max depth 128), the changes on this branch collectively reduce the time required for layout creation by ~50%, and reduce the overall runtime of the end-to-end protocol by ~15%. Your mileage will of course vary, however, and I would anticipate higher relative savings for larger/deeper experiment designs. I have not done profiling for the two-qubit case yet, so I won't posit any guesses on the impact there (I would guess it to be less relatively speaking, but can't easily ballpark this without actually testing).

### Summary of Changes:
Here is a high-level summary of the changes. I will also be leaving some comments on particular bits of code to explain any more detailed reasoning for some changes (when in doubt the answer is likely that I tried a few different things and picked the fastest way to do something).

**Note on synthetic benchmarks:** For the quoted figures below these were mostly done using a randomly generated set of test circuits generated according to the following code snippet:

```
import pygsti
from pygsti.circuits.circuitconstruction import list_random_circuits_onelen
from pygsti.modelpacks import smq1Q_XYI
import numpy as np
from numpy.random import default_rng
from itertools import product
lens = [1,2,4,8,16,32,64]
op_labels = tuple(smq1Q_XYI.target_model().operations.keys())
random_circuits = []
for length in lens:
    random_circuits.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 1234))
random_circuits_2 = []
for length in lens:
    random_circuits_2.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 12345))
rng = default_rng(1234)
rng.shuffle(random_circuits)
rng = default_rng(12345)
rng.shuffle(random_circuits_2)
random_pairs = list(product(random_circuits,repeat=2))
```

For equality testing this was run against `random_pairs`, for inclusion testing against `random_circuits` and `random_circuits_2`, and for other tests against `random_circuits` unless stated otherwise.

1. Circuit equality testing: The primary change was adding a length check to allow quicker short circuiting. A somewhat faster codepath was also added for comparisons between static circuits, avoiding recomputation of a circuit's tuple representation. In synthetic testing this gave a ~3X speed up.
2. Inclusion testing:  Nothing specific was added here, but the change to circuit equality testing also sped this up by a slightly larger factor of ~3.4.
3. Hashing: Hashes are now precomputed and cached at the point at which a circuit is made static (this is an accepted practice in python for immutable objects, e.g. strings cache their hashes). Related to this change, at staticization time we precompute the tuple representation of the circuit and store this as well. This is the thing which actually gets hashed. Previously the `tup` property was called whenever this was needed and this property call regenenerated the tuple representation from scratch each time, which is expensive. In synthetic tests this resulted in a ~4.5X speedup, but in the context of the GST profiling described above this was closer to 10X (the computational savings should scale with circuit depth, so this isn't too shocking of a difference). As part of this change I have added static assertions to a few methods previously missing them for safety (we don't want to inadvertently invalidate the cached hash).
4. Copying: Previously copying a circuit went through the standard `__init__` pathway to create a copied circuit instance. With the new implementation we utilize a much faster specialized codepath inspired by the previously existing `_fastinit` one. But with a few more microptimizations to take advantage of some of the additional promises we can rely on in the context of copying (which we can't as easily in `_fastinit`). Speedups here range from 3-15X depending on the input and output circuit types. With the low end corresponding to switching static to editable (or vice versa), and the high end keeping the input type unchanged.

This is the end of the changes where I profiled these using synthetic tests. For the rest I did this in the context of GST optimization using cProfile/snakeviz, but didn't carefully record before/after relative changes.

5. Microoptimization galore: There are a bunch of microoptimizations of the sort that usually don't matter except for when  code is particularly hot, and the circuit code is. I won't give an exhaustive list, but this includes things like:
    - Attribute access is faster than method/property access/resolution (by 2X or more in some cases). So many trivial (and some non-trivial) property calls have had the corresponding attributes inlined.
    - Inlining functions to reduce function call overhead.
    - List comprehensions are faster than generator expressions. So, e.g. `tuple([i for i in iterable])` is measurably faster than `tuple(i for i in iterable)` in practice.
    - List comprehensions are faster than the 'build via appends in for loop' approach. Often by a notable amount.
    - Using falsy/truthiness directly is faster than the `if len(obj)>0` idiom if all you're doing is checking if an iterable is empty.

6. Rework SeparatePOVMCircuit internals: The `full_effect_labels` property was getting called a lot, and each time it was it was regenerating the result from scratch (with some non-trivial preprocessing). I've updated the internals to cache this, and also deal with updating the cached value when necessary. This gave a measurable speedup.
7. New `sandwich` method for circuits. This takes as input two tuples of Labels which are prepended and appended to a circuit. This is around twice as fast as doing this with a pair of calls to `__add__` (since we only need 1 call to `_fastinit`). This is not an uncommon task (I added as part of some changes to the circuit completion code in OpModel), so if you spot any places where this is being done it'd be worth switching (maybe @rileyjmurray can write a fancy regex to find these instances).

Additional changes (either outside of the Circuit class, or otherwise not directly performance focused):

8. Bulk versions of `complete_circuit` and `split_circuit`: There are now versions of the `complete_circuit` and `split_circuit` methods for the OpModel class, which respectively add prep and povm labels to a circuit or split them off, which operate on a complete list of circuits all at once. This allows for avoiding some model specific lookups (for values which don't change circuit-to-circuit) many times and is faster overall.
9. `expand_subcircuits_inplace`: The implementation of `expand_subcircuits_inplace` has been reworked so that the behavior is now to recursively expand CircuitLabel objects. Previously the expansion routine only took a single pass through the circuit, meaning that if the circuit contained any nested CircuitLabels (which are perfectly happy with being nested, an easy way to get something like this is to call `to_label` on a circuit containing CircuitLabels) we would only end up expanding out the top level. I noticed there were a few places in the layout creation/forward simulation code that implicitly assumed only a single-level of CircuitLabels would be present, which would likely have eventually led to some really annoying bugs to track down. This change covers that edge case, and should only be slightly more expensive than the previous routine in the typical case where there is only a single-level of CircuitLabels.
10. I've added a deprecation warning to the `simulate` method. If there is strong opposition to removing this functionality we can discuss this, of course. My reasoning here is that given the rest of pygsti's structure this functionality really isn't something we should be asking a Circuit to do, as the preferred pattern nowadays is delegate this job to the Model/ForwardSimulator, rather than asking the circuit itself. I don't know how widespread the use of this is (my guess is that this is largely vestigial), but I'd recommend slating this for removal in 0.9.14. Related note, I don't know if we have a running list of things we've deprecated and their expiry date, but if not we should assemble one.

Finally, there are a number of miscellaneous documentation fixes, bugfixes, and updates to unit tests.

P.S. There are some more layout creation performance enhancements inbound on a separate branch, but I've opted to keep these changes together in a separate PR to keep things thematic and make it easier to review without context switching.
Faster Circuit Primitives
I have always been slightly perturbed by the amount of time that things such as COPA layout creation take when performing GST analysis. There are a number of reasons for this, but one family of reasons is that layout generation relies heavily on a number of Circuit primitives that are not as performant as they can be. This PR is a step toward rectifying this through re-implementation of a number of methods identified as bottlenecks through profiling in the context of actual end-to-end GST analysis. 

In terms of the top-level numbers, for 1-qubit GST fitting using the 'full TP' parameterization with an experiment design containing 1120 circuits (max depth 128), the changes on this branch collectively reduce the time required for layout creation by ~50%, and reduce the overall runtime of the end-to-end protocol by ~15%. Your mileage will of course vary, however, and I would anticipate higher relative savings for larger/deeper experiment designs. I have not done profiling for the two-qubit case yet, so I won't posit any guesses on the impact there (I would guess it to be less relatively speaking, but can't easily ballpark this without actually testing).

### Summary of Changes:
Here is a high-level summary of the changes. I will also be leaving some comments on particular bits of code to explain any more detailed reasoning for some changes (when in doubt the answer is likely that I tried a few different things and picked the fastest way to do something).

**Note on synthetic benchmarks:** For the quoted figures below these were mostly done using a randomly generated set of test circuits generated according to the following code snippet:

```
import pygsti
from pygsti.circuits.circuitconstruction import list_random_circuits_onelen
from pygsti.modelpacks import smq1Q_XYI
import numpy as np
from numpy.random import default_rng
from itertools import product
lens = [1,2,4,8,16,32,64]
op_labels = tuple(smq1Q_XYI.target_model().operations.keys())
random_circuits = []
for length in lens:
    random_circuits.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 1234))
random_circuits_2 = []
for length in lens:
    random_circuits_2.extend(list_random_circuits_onelen(op_labels, length, 100, seed = 12345))
rng = default_rng(1234)
rng.shuffle(random_circuits)
rng = default_rng(12345)
rng.shuffle(random_circuits_2)
random_pairs = list(product(random_circuits,repeat=2))
```

For equality testing this was run against `random_pairs`, for inclusion testing against `random_circuits` and `random_circuits_2`, and for other tests against `random_circuits` unless stated otherwise.

1. Circuit equality testing: The primary change was adding a length check to allow quicker short circuiting. A somewhat faster codepath was also added for comparisons between static circuits, avoiding recomputation of a circuit's tuple representation. In synthetic testing this gave a ~3X speed up.
2. Inclusion testing:  Nothing specific was added here, but the change to circuit equality testing also sped this up by a slightly larger factor of ~3.4.
3. Hashing: Hashes are now precomputed and cached at the point at which a circuit is made static (this is an accepted practice in python for immutable objects, e.g. strings cache their hashes). Related to this change, at staticization time we precompute the tuple representation of the circuit and store this as well. This is the thing which actually gets hashed. Previously the `tup` property was called whenever this was needed and this property call regenenerated the tuple representation from scratch each time, which is expensive. In synthetic tests this resulted in a ~4.5X speedup, but in the context of the GST profiling described above this was closer to 10X (the computational savings should scale with circuit depth, so this isn't too shocking of a difference). As part of this change I have added static assertions to a few methods previously missing them for safety (we don't want to inadvertently invalidate the cached hash).
4. Copying: Previously copying a circuit went through the standard `__init__` pathway to create a copied circuit instance. With the new implementation we utilize a much faster specialized codepath inspired by the previously existing `_fastinit` one. But with a few more microptimizations to take advantage of some of the additional promises we can rely on in the context of copying (which we can't as easily in `_fastinit`). Speedups here range from 3-15X depending on the input and output circuit types. With the low end corresponding to switching static to editable (or vice versa), and the high end keeping the input type unchanged.

This is the end of the changes where I profiled these using synthetic tests. For the rest I did this in the context of GST optimization using cProfile/snakeviz, but didn't carefully record before/after relative changes.

5. Microoptimization galore: There are a bunch of microoptimizations of the sort that usually don't matter except for when  code is particularly hot, and the circuit code is. I won't give an exhaustive list, but this includes things like:
    - Attribute access is faster than method/property access/resolution (by 2X or more in some cases). So many trivial (and some non-trivial) property calls have had the corresponding attributes inlined.
    - Inlining functions to reduce function call overhead.
    - List comprehensions are faster than generator expressions. So, e.g. `tuple([i for i in iterable])` is measurably faster than `tuple(i for i in iterable)` in practice.
    - List comprehensions are faster than the 'build via appends in for loop' approach. Often by a notable amount.
    - Using falsy/truthiness directly is faster than the `if len(obj)>0` idiom if all you're doing is checking if an iterable is empty.

6. Rework SeparatePOVMCircuit internals: The `full_effect_labels` property was getting called a lot, and each time it was it was regenerating the result from scratch (with some non-trivial preprocessing). I've updated the internals to cache this, and also deal with updating the cached value when necessary. This gave a measurable speedup.
7. New `sandwich` method for circuits. This takes as input two tuples of Labels which are prepended and appended to a circuit. This is around twice as fast as doing this with a pair of calls to `__add__` (since we only need 1 call to `_fastinit`). This is not an uncommon task (I added as part of some changes to the circuit completion code in OpModel), so if you spot any places where this is being done it'd be worth switching (maybe @rileyjmurray can write a fancy regex to find these instances).

Additional changes (either outside of the Circuit class, or otherwise not directly performance focused):

8. Bulk versions of `complete_circuit` and `split_circuit`: There are now versions of the `complete_circuit` and `split_circuit` methods for the OpModel class, which respectively add prep and povm labels to a circuit or split them off, which operate on a complete list of circuits all at once. This allows for avoiding some model specific lookups (for values which don't change circuit-to-circuit) many times and is faster overall.
9. `expand_subcircuits_inplace`: The implementation of `expand_subcircuits_inplace` has been reworked so that the behavior is now to recursively expand CircuitLabel objects. Previously the expansion routine only took a single pass through the circuit, meaning that if the circuit contained any nested CircuitLabels (which are perfectly happy with being nested, an easy way to get something like this is to call `to_label` on a circuit containing CircuitLabels) we would only end up expanding out the top level. I noticed there were a few places in the layout creation/forward simulation code that implicitly assumed only a single-level of CircuitLabels would be present, which would likely have eventually led to some really annoying bugs to track down. This change covers that edge case, and should only be slightly more expensive than the previous routine in the typical case where there is only a single-level of CircuitLabels.
10. I've added a deprecation warning to the `simulate` method. If there is strong opposition to removing this functionality we can discuss this, of course. My reasoning here is that given the rest of pygsti's structure this functionality really isn't something we should be asking a Circuit to do, as the preferred pattern nowadays is delegate this job to the Model/ForwardSimulator, rather than asking the circuit itself. I don't know how widespread the use of this is (my guess is that this is largely vestigial), but I'd recommend slating this for removal in 0.9.14. Related note, I don't know if we have a running list of things we've deprecated and their expiry date, but if not we should assemble one.

Finally, there are a number of miscellaneous documentation fixes, bugfixes, and updates to unit tests.

P.S. There are some more layout creation performance enhancements inbound on a separate branch, but I've opted to keep these changes together in a separate PR to keep things thematic and make it easier to review without context switching.


Discussion / hints:

Fix patch:
diff --git a/pygsti/algorithms/core.py b/pygsti/algorithms/core.py
index f2b749136..dd0a21ef7 100644
--- a/pygsti/algorithms/core.py
+++ b/pygsti/algorithms/core.py
@@ -33,6 +33,8 @@
 from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
 from pygsti.optimize.customlm import CustomLMOptimizer as _CustomLMOptimizer
 from pygsti.optimize.customlm import Optimizer as _Optimizer
+from pygsti import forwardsims as _fwdsims
+from pygsti import layouts as _layouts
 
 _dummy_profiler = _DummyProfiler()
 
@@ -400,7 +402,7 @@ def _construct_ab(prep_fiducials, effect_fiducials, model, dataset, op_label_ali
         for j, rhostr in enumerate(prep_fiducials):
             opLabelString = rhostr + estr  # LEXICOGRAPHICAL VS MATRIX ORDER
             dsStr = opLabelString.replace_layers_with_aliases(op_label_aliases)
-            expd_circuit_outcomes = opLabelString.expand_instruments_and_separate_povm(model)
+            expd_circuit_outcomes = model.expand_instruments_and_separate_povm(opLabelString)
             assert(len(expd_circuit_outcomes) == 1), "No instruments are allowed in LGST fiducials!"
             unique_key = next(iter(expd_circuit_outcomes.keys()))
             outcomes = expd_circuit_outcomes[unique_key]
@@ -429,7 +431,7 @@ def _construct_x_matrix(prep_fiducials, effect_fiducials, model, op_label_tuple,
         for j, rhostr in enumerate(prep_fiducials):
             opLabelString = rhostr + _circuits.Circuit(op_label_tuple, line_labels=rhostr.line_labels) + estr
             dsStr = opLabelString.replace_layers_with_aliases(op_label_aliases)
-            expd_circuit_outcomes = opLabelString.expand_instruments_and_separate_povm(model)
+            expd_circuit_outcomes = model.expand_instruments_and_separate_povm(opLabelString)
             dsRow_fractions = dataset[dsStr].fractions
             assert(len(expd_circuit_outcomes) == nVariants)
 
@@ -677,16 +679,10 @@ def run_gst_fit(mdc_store, optimizer, objective_function_builder, verbosity=0):
         if _np.linalg.norm(mdc_store.model.to_vector() - v_cmp) > 1e-6:
             raise ValueError("MPI ERROR: *different* MC2GST start models"
                              " given to different processors!")                   # pragma: no cover
-
-    #MEM from ..baseobjs.profiler import Profiler
-    #MEM debug_prof = Profiler(comm)
-    #MEM debug_prof.print_memory("run_gst_fit1", True)
-
+        
     if objective_function_builder is not None:
         objective_function_builder = _objfns.ObjectiveFunctionBuilder.cast(objective_function_builder)
-        #MEM debug_prof.print_memory("run_gst_fit2", True)
         objective = objective_function_builder.build_from_store(mdc_store, printer)  # (objective is *also* a store)
-        #MEM debug_prof.print_memory("run_gst_fit3", True)
     else:
         assert(isinstance(mdc_store, _objfns.ObjectiveFunction)), \
             "When `objective_function_builder` is None, `mdc_store` must be an objective fn!"
@@ -705,14 +701,8 @@ def run_gst_fit(mdc_store, optimizer, objective_function_builder, verbosity=0):
 
     printer.log("Completed in %.1fs" % (_time.time() - tStart), 1)
 
-    #if target_model is not None:
-    #  target_vec = target_model.to_vector()
-    #  targetErrVec = _objective_func(target_vec)
-    #  return minErrVec, soln_gs, targetErrVec
     profiler.add_time("do_mc2gst: total time", tStart)
-    #TODO: evTree.permute_computation_to_original(minErrVec) #Doesn't work b/c minErrVec is flattened
-    # but maybe best to just remove minErrVec from return value since this isn't very useful
-    # anyway?
+
     return opt_result, objective
 
 
@@ -888,10 +878,30 @@ def _max_array_types(artypes_list):  # get the maximum number of each array type
     #The ModelDatasetCircuitsStore
     printer.log('Precomputing CircuitOutcomeProbabilityArray layouts for each iteration.', 2)
     precomp_layouts = []
+
+    #pre-compute a dictionary caching completed circuits for layout construction performance.
+    unique_circuits = list({ckt for circuit_list in circuit_lists for ckt in circuit_list})
+    if isinstance(mdl.sim, (_fwdsims.MatrixForwardSimulator, _fwdsims.MapForwardSimulator)):
+        precomp_layout_circuit_cache = mdl.sim.create_copa_layout_circuit_cache(unique_circuits, mdl, dataset=dataset)
+    else:
+        precomp_layout_circuit_cache = None
+
     for i, circuit_list in enumerate(circuit_lists):
         printer.log(f'Layout for iteration {i}', 2)
-        precomp_layouts.append(mdl.sim.create_layout(circuit_list, dataset, resource_alloc, array_types, verbosity= printer - 1))
-    
+        precomp_layouts.append(mdl.sim.create_layout(circuit_list, dataset, resource_alloc, array_types, verbosity= printer - 1,
+                                                    layout_creation_circuit_cache = precomp_layout_circuit_cache))
+        
+    #precompute a cache of possible outcome counts for each circuits to accelerate MDC store creation
+    if isinstance(mdl, _models.model.OpModel):
+        if precomp_layout_circuit_cache is not None: #then grab the split circuits from there.
+            expanded_circuit_outcome_list = mdl.bulk_expand_instruments_and_separate_povm(unique_circuits, 
+                                                                                          completed_circuits= precomp_layout_circuit_cache['completed_circuits'].values())
+        else:
+            expanded_circuit_outcome_list = mdl.bulk_expand_instruments_and_separate_povm(unique_circuits)        
+        outcome_count_by_circuit_cache = {ckt: len(outcome_tup) for ckt,outcome_tup in zip(unique_circuits, expanded_circuit_outcome_list)}
+    else:
+        outcome_count_by_circuit_cache = {ckt: mdl.compute_num_outcomes(ckt) for ckt in unique_circuits}
+
     with printer.progress_logging(1):
         for i in range(starting_index, len(circuit_lists)):
             circuitsToEstimate = circuit_lists[i]
@@ -908,7 +918,8 @@ def _max_array_types(artypes_list):  # get the maximum number of each array type
             mdl.basis = start_model.basis  # set basis in case of CPTP constraints (needed?)
             initial_mdc_store = _objfns.ModelDatasetCircuitsStore(mdl, dataset, circuitsToEstimate, resource_alloc,
                                                                   array_types=array_types, verbosity=printer - 1, 
-                                                                  precomp_layout = precomp_layouts[i])
+                                                                  precomp_layout = precomp_layouts[i],
+                                                                  outcome_count_by_circuit=outcome_count_by_circuit_cache)
             mdc_store = initial_mdc_store
 
             for j, obj_fn_builder in enumerate(iteration_objfn_builders):
diff --git a/pygsti/circuits/__init__.py b/pygsti/circuits/__init__.py
index 46be7652f..28cd30f66 100644
--- a/pygsti/circuits/__init__.py
+++ b/pygsti/circuits/__init__.py
@@ -10,7 +10,7 @@
 # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
 #***************************************************************************************************
 
-from .circuit import Circuit
+from .circuit import Circuit, SeparatePOVMCircuit
 from .circuitlist import CircuitList
 from .circuitstructure import CircuitPlaquette, FiducialPairPlaquette, \
     GermFiducialPairPlaquette, PlaquetteGridCircuitStructure
diff --git a/pygsti/circuits/circuit.py b/pygsti/circuits/circuit.py
index 2977dddac..2ccb7aaae 100644
--- a/pygsti/circuits/circuit.py
+++ b/pygsti/circuits/circuit.py
@@ -1095,13 +1095,12 @@ def _proc_lines_arg(self, lines):
     def _proc_key_arg(self, key):
         """ Pre-process the key argument used by many methods """
         if isinstance(key, tuple):
-            if len(key) != 2: return IndexError("Index must be of the form <layerIndex>,<lineIndex>")
-            layers = key[0]
-            lines = key[1]
+            if len(key) != 2: 
+                return IndexError("Index must be of the form <layerIndex>,<lineIndex>")
+            else:
+                return key[0], key[1]
         else:
-            layers = key
-            lines = None
-        return layers, lines
+            return key, None
 
     def _layer_components(self, ilayer):
         """ Get the components of the `ilayer`-th layer as a list/tuple. """
@@ -1191,22 +1190,41 @@ def extract_labels(self, layers=None, lines=None, strict=True):
             `layers` is a single integer and as a `Circuit` otherwise.
             Note: if you want a `Circuit` when only selecting one layer,
             set `layers` to a slice or tuple containing just a single index.
+            Note that the returned circuit doesn't retain any original
+            metadata, such as the compilable layer indices or occurence id.
         """
         nonint_layers = not isinstance(layers, int)
 
         #Shortcut for common case when lines == None and when we're only taking a layer slice/index
-        if lines is None:
-            assert(layers is not None)
-            if nonint_layers is False: return self.layertup[layers]
-            if isinstance(layers, slice) and strict is True:  # if strict=False, then need to recompute line labels
-                return Circuit._fastinit(self._labels[layers], self._line_labels, not self._static)
+        if lines is None and layers is not None:
+            if self._static:
+                if not nonint_layers:
+                    return self._labels[layers]
+                if isinstance(layers, slice) and strict is True:  # if strict=False, then need to recompute line labels
+                    #can speed this up a measurably by manually computing the new hashable tuple value and hash
+                    if not self._line_labels in (('*',), ()):
+                        new_hashable_tup = self._labels[layers] + ('@',) + self._line_labels
+                    else:
+                        new_hashable_tup = self._labels[layers]
+                    ret = Circuit.__new__(Circuit)
+                    return ret._copy_init(self._labels[layers], self._line_labels, not self._static, hashable_tup= new_hashable_tup, precomp_hash=hash(new_hashable_tup))
+            else:
+                if not nonint_layers:
+                    return self.layertup[layers]
+                if isinstance(layers, slice) and strict is True:  # if strict=False, then need to recompute line labels
+                    return Circuit._fastinit(self._labels[layers], self._line_labels, not self._static)
+        #otherwise assert both are not None:
+
 
         layers = self._proc_layers_arg(layers)
         lines = self._proc_lines_arg(lines)
         if len(layers) == 0 or len(lines) == 0:
-            return Circuit._fastinit(() if self._static else [],
-                                     tuple(lines) if self._static else lines,
-                                     not self._static) if nonint_layers else None  # zero-area region
+            if self._static:
+                return Circuit._fastinit((), tuple(lines), False)  # zero-area region
+            else:
+                return Circuit._fastinit(() if self._static else [],
+                                        tuple(lines) if self._static else lines,
+                                        not self._static)  # zero-area region
 
         ret = []
         if self._static:
@@ -3447,7 +3465,7 @@ def num_gates(self):
         """
         if self._static:
             def cnt(lbl):  # obj a Label, perhaps compound
-                if lbl.is_simple():  # a simple label
+                if lbl.IS_SIMPLE:  # a simple label
                     return 1 if (lbl.sslbls is not None) else 0
                 else:
                     return sum([cnt(sublbl) for sublbl in lbl.components])
@@ -4427,107 +4445,6 @@ def done_editing(self):
         self._hashable_tup = self.tup
         self._hash = hash(self._hashable_tup)
 
-    def expand_instruments_and_separate_povm(self, model, observed_outcomes=None):
-        """
-        Creates a dictionary of :class:`SeparatePOVMCircuit` objects from expanding the instruments of this circuit.
-
-        Each key of the returned dictionary replaces the instruments in this circuit with a selection
-        of their members.  (The size of the resulting dictionary is the product of the sizes of
-        each instrument appearing in this circuit when `observed_outcomes is None`).  Keys are stored
-        as :class:`SeparatePOVMCircuit` objects so it's easy to keep track of which POVM outcomes (effects)
-        correspond to observed data.  This function is, for the most part, used internally to process
-        a circuit before computing its outcome probabilities.
-
-        Parameters
-        ----------
-        model : Model
-            The model used to provide necessary details regarding the expansion, including:
-
-            - default SPAM layers
-            - definitions of instrument-containing layers
-            - expansions of individual instruments and POVMs
-
-        Returns
-        -------
-        OrderedDict
-            A dict whose keys are :class:`SeparatePOVMCircuit` objects and whose
-            values are tuples of the outcome labels corresponding to this circuit,
-            one per POVM effect held in the key.
-        """
-        complete_circuit = model.complete_circuit(self)
-        expanded_circuit_outcomes = _collections.OrderedDict()
-        povm_lbl = complete_circuit[-1]  # "complete" circuits always end with a POVM label
-        circuit_without_povm = complete_circuit[0:len(complete_circuit) - 1]
-
-        def create_tree(lst):
-            subs = _collections.OrderedDict()
-            for el in lst:
-                if len(el) > 0:
-                    if el[0] not in subs: subs[el[0]] = []
-                    subs[el[0]].append(el[1:])
-            return _collections.OrderedDict([(k, create_tree(sub_lst)) for k, sub_lst in subs.items()])
-
-        def add_expanded_circuit_outcomes(circuit, running_outcomes, ootree, start):
-            """
-            """
-            cir = circuit if start == 0 else circuit[start:]  # for performance, avoid uneeded slicing
-            for k, layer_label in enumerate(cir, start=start):
-                components = layer_label.components
-                #instrument_inds = _np.nonzero([model._is_primitive_instrument_layer_lbl(component)
-                #                               for component in components])[0]  # SLOWER than statement below
-                instrument_inds = _np.array([i for i, component in enumerate(components)
-                                             if model._is_primitive_instrument_layer_lbl(component)])
-                if instrument_inds.size > 0:
-                    # This layer contains at least one instrument => recurse with instrument(s) replaced with
-                    #  all combinations of their members.
-                    component_lookup = {i: comp for i, comp in enumerate(components)}
-                    instrument_members = [model._member_labels_for_instrument(components[i])
-                                          for i in instrument_inds]  # also components of outcome labels
-                    for selected_instrmt_members in _itertools.product(*instrument_members):
-                        expanded_layer_lbl = component_lookup.copy()
-                        expanded_layer_lbl.update({i: components[i] + "_" + sel
-                                                   for i, sel in zip(instrument_inds, selected_instrmt_members)})
-                        expanded_layer_lbl = _Label([expanded_layer_lbl[i] for i in range(len(components))])
-
-                        if ootree is not None:
-                            new_ootree = ootree
-                            for sel in selected_instrmt_members:
-                                new_ootree = new_ootree.get(sel, {})
-                            if len(new_ootree) == 0: continue  # no observed outcomes along this outcome-tree path
-                        else:
-                            new_ootree = None
-
-                        add_expanded_circuit_outcomes(circuit[0:k] + Circuit((expanded_layer_lbl,)) + circuit[k + 1:],
-                                                      running_outcomes + selected_instrmt_members, new_ootree, k + 1)
-                    break
-
-            else:  # no more instruments to process: `cir` contains no instruments => add an expanded circuit
-                assert(circuit not in expanded_circuit_outcomes)  # shouldn't be possible to generate duplicates...
-                elabels = model._effect_labels_for_povm(povm_lbl) if (observed_outcomes is None) \
-                    else tuple(ootree.keys())
-                outcomes = tuple((running_outcomes + (elabel,) for elabel in elabels))
-                expanded_circuit_outcomes[SeparatePOVMCircuit(circuit, povm_lbl, elabels)] = outcomes
-
-        ootree = create_tree(observed_outcomes) if observed_outcomes is not None else None  # tree of observed outcomes
-        # e.g. [('0','00'), ('0','01'), ('1','10')] ==> {'0': {'00': {}, '01': {}}, '1': {'10': {}}}
-
-        if model._has_instruments():
-            add_expanded_circuit_outcomes(circuit_without_povm, (), ootree, start=0)
-        else:
-            # It may be helpful to cache the set of elabels for a POVM (maybe within the model?) because
-            # currently the call to _effect_labels_for_povm may be a bottleneck.  It's needed, even when we have
-            # observed outcomes, because there may be some observed outcomes that aren't modeled (e.g. leakage states)
-            if observed_outcomes is None:
-                elabels = model._effect_labels_for_povm(povm_lbl)
-            else:
-                possible_lbls = set(model._effect_labels_for_povm(povm_lbl))
-                elabels = tuple([oo for oo in ootree.keys() if oo in possible_lbls])
-            outcomes = tuple(((elabel,) for elabel in elabels))
-            expanded_circuit_outcomes[SeparatePOVMCircuit(circuit_without_povm, povm_lbl, elabels)] = outcomes
-
-        return expanded_circuit_outcomes
-
-
 class CompressedCircuit(object):
     """
     A "compressed" Circuit that requires less disk space.
diff --git a/pygsti/data/dataset.py b/pygsti/data/dataset.py
index 6214fbf41..ce7bb52c6 100644
--- a/pygsti/data/dataset.py
+++ b/pygsti/data/dataset.py
@@ -296,34 +296,10 @@ def timeseries_for_outcomes(self):
         last_time = None
         seriesDict = {self.dataset.olIndex[ol]: [] for ol in self.dataset.outcome_labels}
 
-        #REMOVED: (though this gives slightly different behavior)
-        #for outcome_label in self.outcomes:
-        #    if outcome_label not in seriesDict.keys():
-        #        seriesDict[outcome_label] = []
-
         if self.reps is None:
             reps = _np.ones(len(self.time), _np.int64)
         else: reps = self.reps
 
-        # An alternate implementation that appears to be (surprisingly?) slower...
-        ##Get time bin locations
-        #time_bins_borders = []
-        #last_time = None
-        #for i, t in enumerate(self.time):
-        #    if t != last_time:
-        #        time_bins_borders.append(i)
-        #        last_time = t
-        #time_bins_borders.append(len(self.time))
-        #nTimes = len(time_bins_borders) - 1
-        #
-        #seriesDict = {self.dataset.olIndex[ol]: _np.zeros(nTimes, _np.int64) for ol in self.dataset.outcome_labels}
-        #
-        #for i in range(nTimes):
-        #    slc = slice(time_bins_borders[i],time_bins_borders[i+1])
-        #    times.append( self.time[slc.start] )
-        #    for oli, rep in zip(self.oli[slc], reps[slc]):
-        #        seriesDict[oli][i] += rep
-
         for t, oli, rep in zip(self.time, self.oli, reps):
 
             if t != last_time:
@@ -586,26 +562,28 @@ def _get_counts(self, timestamp=None, all_outcomes=False):
             tslc = _np.where(_np.isclose(self.time, timestamp))[0]
         else: tslc = slice(None)
 
+        oli_tslc = self.oli[tslc]
         nOutcomes = len(self.dataset.olIndex)
-        nIndices = len(self.oli[tslc])
+        nIndices = len(oli_tslc)
+        
         if nOutcomes <= nIndices or all_outcomes:
             if self.reps is None:
                 for ol, i in self.dataset.olIndex.items():
-                    cnt = float(_np.count_nonzero(_np.equal(self.oli[tslc], i)))
-                    if all_outcomes or cnt > 0:
+                    cnt = float(_np.count_nonzero(_np.equal(oli_tslc, i)))
+                    if cnt > 0 or all_outcomes:
                         cntDict.setitem_unsafe(ol, cnt)
             else:
                 for ol, i in self.dataset.olIndex.items():
-                    inds = _np.nonzero(_np.equal(self.oli[tslc], i))[0]
-                    if all_outcomes or len(inds) > 0:
+                    inds = oli_tslc[oli_tslc == i]
+                    if len(inds) > 0 or all_outcomes:
                         cntDict.setitem_unsafe(ol, float(sum(self.reps[tslc][inds])))
         else:
             if self.reps is None:
-                for ol_index in self.oli[tslc]:
+                for ol_index in oli_tslc:
                     ol = self.dataset.ol[ol_index]
                     cntDict.setitem_unsafe(ol, 1.0 + cntDict.getitem_unsafe(ol, 0.0))
             else:
-                for ol_index, reps in zip(self.oli[tslc], self.reps[tslc]):
+                for ol_index, reps in zip(oli_tslc, self.reps[tslc]):
                     ol = self.dataset.ol[ol_index]
                     cntDict.setitem_unsafe(ol, reps + cntDict.getitem_unsafe(ol, 0.0))
 
@@ -616,7 +594,8 @@ def counts(self):
         """
         Dictionary of per-outcome counts.
         """
-        if self._cntcache: return self._cntcache  # if not None *and* len > 0
+        if self._cntcache: 
+            return self._cntcache  # if not None *and* len > 0
         ret = self._get_counts()
         if self._cntcache is not None:  # == and empty dict {}
             self._cntcache.update(ret)
@@ -1199,10 +1178,10 @@ def _get_row(self, circuit):
         circuit = _cir.Circuit.cast(circuit)
 
         #Note: cirIndex value is either an int (non-static) or a slice (static)
-        repData = self.repData[self.cirIndex[circuit]] \
-            if (self.repData is not None) else None
-        return _DataSetRow(self, self.oliData[self.cirIndex[circuit]],
-                           self.timeData[self.cirIndex[circuit]], repData,
+        cirIndex = self.cirIndex[circuit]
+        repData = self.repData[cirIndex] if (self.repData is not None) else None
+        return _DataSetRow(self, self.oliData[cirIndex],
+                           self.timeData[cirIndex], repData,
                            self.cnt_cache[circuit] if self.bStatic else None,
                            self.auxInfo[circuit])
 
diff --git a/pygsti/forwardsims/forwardsim.py b/pygsti/forwardsims/forwardsim.py
index c5e61b057..2ae19f2f3 100644
--- a/pygsti/forwardsims/forwardsim.py
+++ b/pygsti/forwardsims/forwardsim.py
@@ -323,7 +323,8 @@ def hprobs(self, circuit, resource_alloc=None):
     # ---------------------------------------------------------------------------
 
     def create_layout(self, circuits, dataset=None, resource_alloc=None,
-                      array_types=(), derivative_dimensions=None, verbosity=0):
+                      array_types=(), derivative_dimensions=None, verbosity=0,
+                      layout_creation_circuit_cache = None):
         """
         Constructs an circuit-outcome-probability-array (COPA) layout for `circuits` and `dataset`.
 
@@ -364,6 +365,8 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None,
         verbosity : int or VerbosityPrinter
             Determines how much output to send to stdout.  0 means no output, higher
             integers mean more output.
+        
+        
 
         Returns
         -------
@@ -378,6 +381,15 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None,
                 derivative_dimensions = tuple()
         return _CircuitOutcomeProbabilityArrayLayout.create_from(circuits, self.model, dataset, derivative_dimensions,
                                                                  resource_alloc=resource_alloc)
+    
+    @staticmethod
+    def create_copa_layout_circuit_cache(circuits, model, dataset=None):
+        """
+        Helper function for pre-computing/pre-processing circuits structures
+        used in matrix layout creation.
+        """
+        msg = "Not currently implemented for this forward simulator class."
+        raise NotImplementedError(msg)
 
     def bulk_probs(self, circuits, clip_to=None, resource_alloc=None, smartc=None):
         """
diff --git a/pygsti/forwardsims/mapforwardsim.py b/pygsti/forwardsims/mapforwardsim.py
index 6b19e8d39..6e5ed4f83 100644
--- a/pygsti/forwardsims/mapforwardsim.py
+++ b/pygsti/forwardsims/mapforwardsim.py
@@ -26,6 +26,8 @@
 from pygsti.tools import sharedmemtools as _smt
 from pygsti.tools import slicetools as _slct
 from pygsti.tools.matrixtools import _fas
+from pygsti.tools import listtools as _lt
+from pygsti.circuits import CircuitList as _CircuitList
 
 _dummy_profiler = _DummyProfiler()
 
@@ -49,7 +51,7 @@ class SimpleMapForwardSimulator(_ForwardSimulator):
     # If this is done, then MapForwardSimulator wouldn't need to separately subclass DistributableForwardSimulator
 
     def _compute_circuit_outcome_probabilities(self, array_to_fill, circuit, outcomes, resource_alloc, time=None):
-        expanded_circuit_outcomes = circuit.expand_instruments_and_separate_povm(self.model, outcomes)
+        expanded_circuit_outcomes = self.model.expand_instruments_and_separate_povm(circuit, outcomes)
         outcome_to_index = {outc: i for i, outc in enumerate(outcomes)}
         for spc, spc_outcomes in expanded_circuit_outcomes.items():  # spc is a SeparatePOVMCircuit
             # Note: `spc.circuit_without_povm` *always* begins with a prep label.
@@ -193,7 +195,7 @@ def copy(self):
                                    self._processor_grid, self._pblk_sizes)
 
     def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types=('E',),
-                      derivative_dimensions=None, verbosity=0):
+                      derivative_dimensions=None, verbosity=0, layout_creation_circuit_cache=None):
         """
         Constructs an circuit-outcome-probability-array (COPA) layout for a list of circuits.
 
@@ -223,6 +225,11 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
         verbosity : int or VerbosityPrinter
             Determines how much output to send to stdout.  0 means no output, higher
             integers mean more output.
+        
+        A precomputed dictionary serving as a cache for completed
+            circuits. I.e. circuits with prep labels and POVM labels appended.
+            Along with other useful pre-computed circuit structures used in layout
+            creation.
 
         Returns
         -------
@@ -265,7 +272,8 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
         assert(_np.prod((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
 
         layout = _MapCOPALayout(circuits, self.model, dataset, self._max_cache_size, natoms, na, npp,
-                                param_dimensions, param_blk_sizes, resource_alloc, verbosity)
+                                param_dimensions, param_blk_sizes, resource_alloc, verbosity, 
+                                layout_creation_circuit_cache= layout_creation_circuit_cache)
 
         if mem_limit is not None:
             loc_nparams1 = num_params / npp[0] if len(npp) > 0 else 0
@@ -306,6 +314,41 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
                 printer.log("   Esimated memory required = %.1fGB" % (mem_estimate * GB))
 
         return layout
+    
+    @staticmethod
+    def create_copa_layout_circuit_cache(circuits, model, dataset=None):
+        """
+        Helper function for pre-computing/pre-processing circuits structures
+        used in matrix layout creation.
+        """
+        cache = dict()
+        completed_circuits = model.complete_circuits(circuits)
+
+        cache['completed_circuits'] = {ckt: comp_ckt for ckt, comp_ckt in zip(circuits, completed_circuits)}
+
+        split_circuits = model.split_circuits(completed_circuits, split_prep=False)    
+        cache['split_circuits'] = {ckt: split_ckt for ckt, split_ckt in zip(circuits, split_circuits)}
+        
+
+        if dataset is not None:
+            aliases = circuits.op_label_aliases if isinstance(circuits, _CircuitList) else None
+            ds_circuits = _lt.apply_aliases_to_circuits(circuits, aliases)
+            unique_outcomes_list = []
+            for ckt in ds_circuits:
+                ds_row = dataset[ckt]
+                unique_outcomes_list.append(ds_row.unique_outcomes if ds_row is not None else None)
+        else:
+            unique_outcomes_list = [None]*len(circuits)
+
+        expanded_circuit_outcome_list = model.bulk_expand_instruments_and_separate_povm(circuits, 
+                                                                                        observed_outcomes_list = unique_outcomes_list, 
+                                                                                        completed_circuits= completed_circuits)
+        
+        expanded_circuit_cache = {ckt: expanded_ckt for ckt,expanded_ckt in zip(completed_circuits, expanded_circuit_outcome_list)}                
+        cache['expanded_and_separated_circuits'] = expanded_circuit_cache
+
+        return cache
+
 
     def _bulk_fill_probs_atom(self, array_to_fill, layout_atom, resource_alloc):
         # Note: *don't* set dest_indices arg = layout.element_slice, as this is already done by caller
diff --git a/pygsti/forwardsims/matrixforwardsim.py b/pygsti/forwardsims/matrixforwardsim.py
index ddc18270a..a24c1322d 100644
--- a/pygsti/forwardsims/matrixforwardsim.py
+++ b/pygsti/forwardsims/matrixforwardsim.py
@@ -29,6 +29,8 @@
 from pygsti.tools import sharedmemtools as _smt
 from pygsti.tools import slicetools as _slct
 from pygsti.tools.matrixtools import _fas
+from pygsti.tools import listtools as _lt
+from pygsti.circuits import CircuitList as _CircuitList
 
 _dummy_profiler = _DummyProfiler()
 
@@ -590,7 +592,7 @@ def _compute_circuit_outcome_probabilities(self, array_to_fill, circuit, outcome
         use_scaling = False  # Hardcoded for now
         assert(time is None), "MatrixForwardSimulator cannot be used to simulate time-dependent circuits"
 
-        expanded_circuit_outcomes = circuit.expand_instruments_and_separate_povm(self.model, outcomes)
+        expanded_circuit_outcomes = self.model.expand_instruments_and_separate_povm(circuit, outcomes)
         outcome_to_index = {outc: i for i, outc in enumerate(outcomes)}
         for spc, spc_outcomes in expanded_circuit_outcomes.items():  # spc is a SeparatePOVMCircuit
             indices = [outcome_to_index[o] for o in spc_outcomes]
@@ -1025,7 +1027,7 @@ def _compute_hproduct_cache(self, layout_atom_tree, prod_cache, d_prod_cache1,
         return hProdCache
 
     def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types=('E',),
-                      derivative_dimensions=None, verbosity=0):
+                      derivative_dimensions=None, verbosity=0, layout_creation_circuit_cache= None):
         """
         Constructs an circuit-outcome-probability-array (COPA) layout for a list of circuits.
 
@@ -1056,6 +1058,12 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
             Determines how much output to send to stdout.  0 means no output, higher
             integers mean more output.
 
+        layout_creation_circuit_cache : dict, optional (default None)
+            A precomputed dictionary serving as a cache for completed
+            circuits. I.e. circuits with prep labels and POVM labels appended.
+            Along with other useful pre-computed circuit structures used in layout
+            creation.
+            
         Returns
         -------
         MatrixCOPALayout
@@ -1105,7 +1113,8 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
         assert(_np.prod((na,) + npp) <= nprocs), "Processor grid size exceeds available processors!"
 
         layout = _MatrixCOPALayout(circuits, self.model, dataset, natoms,
-                                   na, npp, param_dimensions, param_blk_sizes, resource_alloc, verbosity)
+                                   na, npp, param_dimensions, param_blk_sizes, resource_alloc, verbosity, 
+                                   layout_creation_circuit_cache=layout_creation_circuit_cache)
 
         if mem_limit is not None:
             loc_nparams1 = num_params / npp[0] if len(npp) > 0 else 0
@@ -1147,6 +1156,46 @@ def create_layout(self, circuits, dataset=None, resource_alloc=None, array_types
                 printer.log("   Esimated memory required = %.1fGB" % (mem_estimate * GB))
 
         return layout
+    
+    @staticmethod
+    def create_copa_layout_circuit_cache(circuits, model, dataset=None):
+        """
+        Helper function for pre-computing/pre-processing circuits structures
+        used in matrix layout creation.
+        """
+        cache = dict()
+        completed_circuits, split_circuits = model.complete_circuits(circuits, return_split=True)
+
+        cache['completed_circuits'] = {ckt: comp_ckt for ckt, comp_ckt in zip(circuits, completed_circuits)}
+        cache['split_circuits'] = {ckt: split_ckt for ckt, split_ckt in zip(circuits, split_circuits)}
+
+        if dataset is not None:
+            aliases = circuits.op_label_aliases if isinstance(circuits, _CircuitList) else None
+            ds_circuits = _lt.apply_aliases_to_circuits(circuits, aliases)
+            unique_outcomes_list = []
+            for ckt in ds_circuits:
+                ds_row = dataset[ckt]
+                unique_outcomes_list.append(ds_row.unique_outcomes if ds_row is not None else None)
+        else:
+            unique_outcomes_list = [None]*len(circuits)
+
+        expanded_circuit_outcome_list = model.bulk_expand_instruments_and_separate_povm(circuits, 
+                                                                                        observed_outcomes_list = unique_outcomes_list, 
+                                                                                        split_circuits = split_circuits)
+        
+        expanded_circuit_cache = {ckt: expanded_ckt for ckt,expanded_ckt in zip(circuits, expanded_circuit_outcome_list)}
+                    
+        cache['expanded_and_separated_circuits'] = expanded_circuit_cache
+
+        expanded_subcircuits_no_spam_cache = dict()
+        for expc_outcomes in cache['expanded_and_separated_circuits'].values():
+            for sep_povm_c, _ in expc_outcomes.items():  # for each expanded cir from unique_i-th circuit
+                exp_nospam_c = sep_povm_c.circuit_without_povm[1:] 
+                expanded_subcircuits_no_spam_cache[exp_nospam_c] = exp_nospam_c.expand_subcircuits()
+
+        cache['expanded_subcircuits_no_spam'] = expanded_subcircuits_no_spam_cache
+
+        return cache
 
     def _scale_exp(self, scale_exps):
         old_err = _np.seterr(over='ignore')
diff --git a/pygsti/forwardsims/torchfwdsim.py b/pygsti/forwardsims/torchfwdsim.py
index 1285e51de..5d61c3f0d 100644
--- a/pygsti/forwardsims/torchfwdsim.py
+++ b/pygsti/forwardsims/torchfwdsim.py
@@ -74,8 +74,9 @@ class StatelessModel:
     def __init__(self, model: ExplicitOpModel, layout: CircuitOutcomeProbabilityArrayLayout):
         circuits = []
         self.outcome_probs_dim = 0
+        #TODO: Refactor this to use the bulk_expand_instruments_and_separate_povm codepath
         for _, circuit, outcomes in layout.iter_unique_circuits():
-            expanded_circuits = circuit.expand_instruments_and_separate_povm(model, outcomes)
+            expanded_circuits = model.expand_instruments_and_separate_povm(circuit, outcomes)
             if len(expanded_circuits) > 1:
                 raise NotImplementedError("I don't know what to do with this.")
             spc = next(iter(expanded_circuits))
diff --git a/pygsti/forwardsims/weakforwardsim.py b/pygsti/forwardsims/weakforwardsim.py
index 017973a1e..32d0e4bc6 100644
--- a/pygsti/forwardsims/weakforwardsim.py
+++ b/pygsti/forwardsims/weakforwardsim.py
@@ -55,7 +55,7 @@ def _compute_circuit_outcome_for_shot(self, circuit, resource_alloc, time=None,
         circuit : Circuit
             A tuple-like object of *simplified* gates (e.g. may include
             instrument elements like 'Imyinst_0') generated by
-            Circuit.expand_instruments_and_separate_povm()
+            OpModel.expand_instruments_and_separate_povm()
 
         resource_alloc: ResourceAlloc
             Currently not used
@@ -77,7 +77,7 @@ def _compute_circuit_outcome_for_shot(self, circuit, resource_alloc, time=None,
         assert(resource_alloc is None), "WeakForwardSimulator cannot use a resource_alloc for one shot."
 
         #prep_label, op_labels, povm_label = self.model.split_circuit(spc_circuit)
-        spc_dict = circuit.expand_instruments_and_separate_povm(self.model,
+        spc_dict = self.model.expand_instruments_and_separate_povm(circuit,
                                                                 observed_outcomes=None)  # FUTURE: observed outcomes?
         assert(len(spc_dict) == 1), "Circuits with instruments are not supported by weak forward simulator (yet)"
         spc = next(iter(spc_dict.keys()))  # first & only SeparatePOVMCircuit
diff --git a/pygsti/layouts/copalayout.py b/pygsti/layouts/copalayout.py
index 430fb6734..bd5020aa8 100644
--- a/pygsti/layouts/copalayout.py
+++ b/pygsti/layouts/copalayout.py
@@ -190,24 +190,26 @@ def __init__(self, circuits, unique_circuits, to_unique, elindex_outcome_tuples,
         if unique_circuits is None and to_unique is None:
             unique_circuits, to_unique = self._compute_unique_circuits(circuits)
         self._unique_circuits = unique_circuits
-        self._unique_circuit_index = _collections.OrderedDict(
-            [(c, i) for i, c in enumerate(self._unique_circuits)])  # original circuits => unique circuit indices
+        self._unique_circuit_index = {c:i for i, c in enumerate(self._unique_circuits)}  # original circuits => unique circuit indices
         self._to_unique = to_unique  # original indices => unique circuit indices
         self._unique_complete_circuits = unique_complete_circuits  # Note: can be None
         self._param_dimensions = param_dimensions
         self._resource_alloc = _ResourceAllocation.cast(resource_alloc)
 
-        max_element_index = max(_it.chain(*[[ei for ei, _ in pairs] for pairs in elindex_outcome_tuples.values()])) \
-            if len(elindex_outcome_tuples) > 0 else -1  # -1 makes _size = 0 below
-        indices = set(i for tuples in elindex_outcome_tuples.values() for i, o in tuples)
+        indices = [i for tuples in elindex_outcome_tuples.values() for i, _ in tuples]
+        max_element_index = max(indices) if len(elindex_outcome_tuples) > 0 else -1  # -1 makes _size = 0 below
+        indices = set(indices)
+        
+        
         self._size = max_element_index + 1
         assert(len(indices) == self._size), \
             "Inconsistency: %d distinct indices but max index + 1 is %d!" % (len(indices), self._size)
 
-        self._outcomes = _collections.OrderedDict()
-        self._element_indices = _collections.OrderedDict()
+        self._outcomes = dict()
+        self._element_indices = dict()
+        sort_idx_func = lambda x: x[0]
         for i_unique, tuples in elindex_outcome_tuples.items():
-            sorted_tuples = sorted(tuples, key=lambda x: x[0])  # sort by element index
+            sorted_tuples = sorted(tuples, key=sort_idx_func)  # sort by element index
             elindices, outcomes = zip(*sorted_tuples)  # sorted by elindex so we make slices whenever possible
             self._outcomes[i_unique] = tuple(outcomes)
             self._element_indices[i_unique] = _slct.list_to_slice(elindices, array_ok=True)
diff --git a/pygsti/layouts/distlayout.py b/pygsti/layouts/distlayout.py
index 7a7184529..9db1150d8 100644
--- a/pygsti/layouts/distlayout.py
+++ b/pygsti/layouts/distlayout.py
@@ -360,9 +360,6 @@ def __init__(self, circuits, unique_circuits, to_unique, unique_complete_circuit
             to_send = 0  # default = contribute nothing to MPI.SUM below
 
             if i in atoms_dict:
-                #print("DB (%d): updating elindex_outcome_tuples w/Atom %d:\n%s"
-                #      % (rank, i, "\n".join(["%d: %s" % (indx, str(tups))
-                #                             for indx, tups in atoms_dict[i].elindex_outcome_tuples.items()])))
                 if start is None:
                     start = stop = offset
                 assert(stop == offset)  # This should be checked by _assert_sequential(myAtomIndices) above
@@ -810,42 +807,6 @@ def __init__(self, circuits, unique_circuits, to_unique, unique_complete_circuit
         super().__init__(local_circuits, local_unique_circuits, local_to_unique, local_elindex_outcome_tuples,
                          local_unique_complete_circuits, param_dimensions, resource_alloc)
 
-        #DEBUG LAYOUT PRINTING
-        #def cnt_str(cnt):
-        #    if cnt is None: return "%4s" % '-'
-        #    return "%4d" % cnt
-        #def slc_str(slc):
-        #    if slc is None: return "%14s" % '--'
-        #    return "%3d->%3d (%3d)" % (slc.start, slc.stop, slc.stop - slc.start) \
-        #        if isinstance(slc, slice) else "%14s" % str(slc)
-        #shm = bool(resource_alloc.host_comm is not None)  # shared mem?
-        #if rank == 0:
-        #    print("%11s %-14s %-14s %-14s   %-14s %-4s %-14s %-4s %-14s %-4s" % (
-        #        '#', 'g-elements', 'g-params', 'g-param2s',
-        #        'h-elements','tot', 'h-params','tot', 'h-params2','tot'),
-        #          flush=True)
-        #resource_alloc.comm.barrier()
-        #for r in range(resource_alloc.comm.size):
-        #    if r == rank:
-        #        my_desc = ("%3d (%2d.%2d)" % (rank, resource_alloc.host_index, resource_alloc.host_comm.rank)) \
-        #                  if shm else ("%11d" % rank)
-        #        print(my_desc, slc_str(self.global_element_slice), slc_str(self.global_param_slice),
-        #              slc_str(self.global_param2_slice), ' ',
-        #              slc_str(self.host_element_slice), cnt_str(self.host_num_elements),
-        #              slc_str(self.host_param_slice), cnt_str(self.host_num_params),
-        #              slc_str(self.host_param2_slice), cnt_str(self.host_num_params2),  flush=True)
-        #    resource_alloc.comm.barrier()
-        #
-        #if rank == 0:
-        #    print("%11s %-14s %-14s %-4s" % ('#', 'g-pfine', 'h-pfine', 'tot'), flush=True)
-        #resource_alloc.comm.barrier()
-        #for r in range(resource_alloc.comm.size):
-        #    if r == rank:
-        #        my_desc = ("%3d (%2d.%2d)" % (rank, resource_alloc.host_index, resource_alloc.host_comm.rank)) \
-        #                  if shm else ("%11d" % rank)
-        #        print(my_desc, slc_str(self.global_param_fine_slice), slc_str(self.host_param_fine_slice),
-        #              cnt_str(self.host_num_params_fine), flush=True)
-        #    resource_alloc.comm.barrier()
 
     @property
     def max_atom_elements(self):
diff --git a/pygsti/layouts/maplayout.py b/pygsti/layouts/maplayout.py
index ca5ca642e..b2142a5e0 100644
--- a/pygsti/layouts/maplayout.py
+++ b/pygsti/layouts/maplayout.py
@@ -51,16 +51,22 @@ class _MapCOPALayoutAtom(_DistributableAtom):
     """
 
     def __init__(self, unique_complete_circuits, ds_circuits, group, model,
-                 dataset, max_cache_size):
+                 dataset, max_cache_size, expanded_complete_circuit_cache = None):
+
+        expanded_circuit_info_by_unique = dict()
+        expanded_circuit_set = dict() # only use SeparatePOVMCircuit keys as ordered set
+
+        if expanded_complete_circuit_cache is None:
+            expanded_complete_circuit_cache = dict()
 
-        expanded_circuit_info_by_unique = _collections.OrderedDict()
-        expanded_circuit_set = _collections.OrderedDict()  # only use SeparatePOVMCircuit keys as ordered set
         for i in group:
-            observed_outcomes = None if (dataset is None) else dataset[ds_circuits[i]].outcomes
-            d = unique_complete_circuits[i].expand_instruments_and_separate_povm(model, observed_outcomes)
+            d = expanded_complete_circuit_cache.get(unique_complete_circuits[i], None)
+            if d is None:
+                unique_observed_outcomes = None if (dataset is None) else dataset[ds_circuits[i]].unique_outcomes
+                d = model.expand_instruments_and_separate_povm(unique_complete_circuits[i], unique_observed_outcomes)
             expanded_circuit_info_by_unique[i] = d  # a dict of SeparatePOVMCircuits => tuples of outcome labels
             expanded_circuit_set.update(d)
-
+            
         expanded_circuits = list(expanded_circuit_set.keys())
         self.table = _PrefixTable(expanded_circuits, max_cache_size)
 
@@ -91,8 +97,7 @@ def __init__(self, unique_complete_circuits, ds_circuits, group, model,
         self.outcomes_by_expcircuit = {}
         self.povm_and_elbls_by_expcircuit = {}
 
-        elindex_outcome_tuples = _collections.OrderedDict([
-            (unique_i, list()) for unique_i in range(len(unique_complete_circuits))])
+        elindex_outcome_tuples = {unique_i: list() for unique_i in range(len(unique_complete_circuits))}
 
         #Assign element indices, "global" indices starting at `offset`
         local_offset = 0
@@ -206,13 +211,42 @@ class MapCOPALayout(_DistributableCOPALayout):
 
     def __init__(self, circuits, model, dataset=None, max_cache_size=None,
                  num_sub_tables=None, num_table_processors=1, num_param_dimension_processors=(),
-                 param_dimensions=(), param_dimension_blk_sizes=(), resource_alloc=None, verbosity=0):
+                 param_dimensions=(), param_dimension_blk_sizes=(), resource_alloc=None, verbosity=0,
+                 layout_creation_circuit_cache=None):
 
         unique_circuits, to_unique = self._compute_unique_circuits(circuits)
         aliases = circuits.op_label_aliases if isinstance(circuits, _CircuitList) else None
         ds_circuits = _lt.apply_aliases_to_circuits(unique_circuits, aliases)
-        unique_complete_circuits = [model.complete_circuit(c) for c in unique_circuits]
-        unique_povmless_circuits = [model.split_circuit(c, split_prep=False)[1] for c in unique_complete_circuits]
+
+        #extract subcaches from layout_creation_circuit_cache:
+        if layout_creation_circuit_cache is None:
+            layout_creation_circuit_cache = dict()
+        self.completed_circuit_cache = layout_creation_circuit_cache.get('completed_circuits', None)
+        self.split_circuit_cache = layout_creation_circuit_cache.get('split_circuits', None)
+        self.expanded_and_separated_circuits_cache = layout_creation_circuit_cache.get('expanded_and_separated_circuits', None)
+        
+        if self.completed_circuit_cache is None:
+            unique_complete_circuits = model.complete_circuits(unique_circuits)
+            split_circuits = model.split_circuits(unique_complete_circuits, split_prep=False)
+        else:
+            unique_complete_circuits = []
+            for c in unique_circuits:
+                comp_ckt = self.completed_circuit_cache.get(c, None)
+                if comp_ckt is not None:
+                    unique_complete_circuits.append(comp_ckt)
+                else:
+                    unique_complete_circuits.append(model.complete_circuit(c))
+            split_circuits = []
+            for c, c_complete in zip(unique_circuits,unique_complete_circuits):
+                split_ckt = self.split_circuit_cache.get(c, None)
+                if split_ckt is not None:
+                    split_circuits.append(split_ckt)
+                else:
+                    split_circuits.append(model.split_circuit(c_complete, split_prep=False))
+            
+
+        #construct list of unique POVM-less circuits.
+        unique_povmless_circuits = [ckt_tup[1] for ckt_tup in split_circuits]
 
         max_sub_table_size = None  # was an argument but never used; remove in future
         if (num_sub_tables is not None and num_sub_tables > 1) or max_sub_table_size is not None:
@@ -221,19 +255,10 @@ def __init__(self, circuits, model, dataset=None, max_cache_size=None,
         else:
             groups = [set(range(len(unique_complete_circuits)))]
 
-        #atoms = []
-        #elindex_outcome_tuples = _collections.OrderedDict(
-        #    [(unique_i, list()) for unique_i in range(len(unique_circuits))])
-
-        #offset = 0
-        #for group in groups:
-        #    atoms.append(_MapCOPALayoutAtom(unique_complete_circuits, ds_circuits, to_orig, group,
-        #                                    model, dataset, offset, elindex_outcome_tuples, max_cache_size))
-        #    offset += atoms[-1].num_elements
-
         def _create_atom(group):
             return _MapCOPALayoutAtom(unique_complete_circuits, ds_circuits, group,
-                                      model, dataset, max_cache_size)
+                                      model, dataset, max_cache_size,
+                                      expanded_complete_circuit_cache=self.expanded_and_separated_circuits_cache)
 
         super().__init__(circuits, unique_circuits, to_unique, unique_complete_circuits,
                          _create_atom, groups, num_table_processors,
@@ -247,4 +272,4 @@ def _create_atom(group):
         unique_to_orig = {unique_i: orig_i for orig_i, unique_i in self._to_unique.items()}  # unique => orig. indices
         for atom in self.atoms:
             for expanded_circuit_i, unique_i in atom.unique_indices_by_expcircuit.items():
-                atom.orig_indices_by_expcircuit[expanded_circuit_i] = unique_to_orig[unique_i]
+                atom.orig_indices_by_expcircuit[expanded_circuit_i] = unique_to_orig[unique_i]
\ No newline at end of file
diff --git a/pygsti/layouts/matrixlayout.py b/pygsti/layouts/matrixlayout.py
index bc92cc054..bfff25a31 100644
--- a/pygsti/layouts/matrixlayout.py
+++ b/pygsti/layouts/matrixlayout.py
@@ -61,6 +61,9 @@ class _MatrixCOPALayoutAtom(_DistributableAtom):
     model : Model
         The model being used to construct this layout.  Used for expanding instruments
         within the circuits.
+    
+    unique_circuits : list of Circuits
+        A list of the unique :class:`Circuit` objects representing the circuits this layout will include.
 
     dataset : DataSet
         The dataset, used to include only observed circuit outcomes in this atom
@@ -68,7 +71,11 @@ class _MatrixCOPALayoutAtom(_DistributableAtom):
     """
 
     def __init__(self, unique_complete_circuits, unique_nospam_circuits, circuits_by_unique_nospam_circuits,
-                 ds_circuits, group, helpful_scratch, model, dataset):
+                 ds_circuits, group, helpful_scratch, model, unique_circuits, dataset=None, expanded_and_separated_circuit_cache=None,
+                 double_expanded_nospam_circuits_cache = None):
+
+        if expanded_and_separated_circuit_cache is None:
+            expanded_and_separated_circuit_cache = dict()
 
         #Note: group gives unique_nospam_circuits indices, which circuits_by_unique_nospam_circuits
         # turns into "unique complete circuit" indices, which the layout via it's to_unique can map
@@ -78,24 +85,24 @@ def add_expanded_circuits(indices, add_to_this_dict):
             for i in indices:
                 nospam_c = unique_nospam_circuits[i]
                 for unique_i in circuits_by_unique_nospam_circuits[nospam_c]:  # "unique" circuits: add SPAM to nospam_c
-                    observed_outcomes = None if (dataset is None) else dataset[ds_circuits[unique_i]].unique_outcomes
-                    expc_outcomes = unique_complete_circuits[unique_i].expand_instruments_and_separate_povm(
-                        model, observed_outcomes)
-                    #Note: unique_complete_circuits may have duplicates (they're only unique *pre*-completion)
-
+                    #the cache is indexed into using the (potentially) incomplete circuits
+                    expc_outcomes = expanded_and_separated_circuit_cache.get(unique_circuits[unique_i], None)
+                    if expc_outcomes is None: #fall back on original non-cache behavior.
+                        observed_outcomes = None if (dataset is None) else dataset[ds_circuits[unique_i]].unique_outcomes
+                        expc_outcomes = model.expand_instruments_and_separate_povm(unique_complete_circuits[unique_i], observed_outcomes)
+                        #and add this new value to the cache.
+                        expanded_and_separated_circuit_cache[unique_circuits[unique_i]] = expc_outcomes 
                     for sep_povm_c, outcomes in expc_outcomes.items():  # for each expanded cir from unique_i-th circuit
                         prep_lbl = sep_povm_c.circuit_without_povm[0]
                         exp_nospam_c = sep_povm_c.circuit_without_povm[1:]  # sep_povm_c *always* has prep lbl
                         spam_tuples = [(prep_lbl, elabel) for elabel in sep_povm_c.full_effect_labels]
-                        outcome_by_spamtuple = _collections.OrderedDict([(st, outcome)
-                                                                         for st, outcome in zip(spam_tuples, outcomes)])
+                        outcome_by_spamtuple = {st:outcome for st, outcome in zip(spam_tuples, outcomes)}
 
                         #Now add these outcomes to `expanded_nospam_circuit_outcomes` - note that multiple "unique_i"'s
                         # may exist for the same expanded & without-spam circuit (exp_nospam_c) and so we need to
                         # keep track of a list of unique_i indices for each circut and spam tuple below.
                         if exp_nospam_c not in _expanded_nospam_circuit_outcomes:
-                            _expanded_nospam_circuit_outcomes[exp_nospam_c] = _collections.OrderedDict(
-                                [(st, (outcome, [unique_i])) for st, outcome in zip(spam_tuples, outcomes)])
+                            _expanded_nospam_circuit_outcomes[exp_nospam_c] = {st:(outcome, [unique_i]) for st, outcome in zip(spam_tuples, outcomes)}
                         else:
                             for st, outcome in outcome_by_spamtuple.items():
                                 if st in _expanded_nospam_circuit_outcomes[exp_nospam_c]:
@@ -109,25 +116,28 @@ def add_expanded_circuits(indices, add_to_this_dict):
 
         # keys = expanded circuits w/out SPAM layers; values = spamtuple => (outcome, unique_is) dictionary that
         # keeps track of which "unique" circuit indices having each spamtuple / outcome.
-        expanded_nospam_circuit_outcomes = _collections.OrderedDict()
+        expanded_nospam_circuit_outcomes = dict()
         add_expanded_circuits(group, expanded_nospam_circuit_outcomes)
-        expanded_nospam_circuits = _collections.OrderedDict(
-            [(i, cir) for i, cir in enumerate(expanded_nospam_circuit_outcomes.keys())])
+        expanded_nospam_circuits = {i:cir for i, cir in enumerate(expanded_nospam_circuit_outcomes.keys())}
 
         # add suggested scratch to the "final" elements as far as the tree creation is concerned
         # - this allows these scratch element to help balance the tree.
-        expanded_nospam_circuit_outcomes_plus_scratch = expanded_nospam_circuit_outcomes.copy()
-        add_expanded_circuits(helpful_scratch, expanded_nospam_circuit_outcomes_plus_scratch)
-        expanded_nospam_circuits_plus_scratch = _collections.OrderedDict(
-            [(i, cir) for i, cir in enumerate(expanded_nospam_circuit_outcomes_plus_scratch.keys())])
-
-        double_expanded_nospam_circuits_plus_scratch = _collections.OrderedDict()
+        if helpful_scratch:
+            expanded_nospam_circuit_outcomes_plus_scratch = expanded_nospam_circuit_outcomes.copy()
+            add_expanded_circuits(helpful_scratch, expanded_nospam_circuit_outcomes_plus_scratch)
+            expanded_nospam_circuits_plus_scratch = {i:cir for i, cir in enumerate(expanded_nospam_circuit_outcomes_plus_scratch.keys())}
+        else:
+            expanded_nospam_circuits_plus_scratch = expanded_nospam_circuits.copy()
+        
+        if double_expanded_nospam_circuits_cache is None:
+            double_expanded_nospam_circuits_cache = dict()
+        double_expanded_nospam_circuits_plus_scratch = dict()
         for i, cir in expanded_nospam_circuits_plus_scratch.items():
-            cir = cir.copy(editable=True)
-            cir.expand_subcircuits()  # expand sub-circuits for a more efficient tree
-            cir.done_editing()
-            double_expanded_nospam_circuits_plus_scratch[i] = cir
-
+            # expand sub-circuits for a more efficient tree
+            double_expanded_ckt = double_expanded_nospam_circuits_cache.get(cir, None)
+            if double_expanded_ckt is None: #Fall back to standard behavior and do expansion.
+                double_expanded_ckt = cir.expand_subcircuits()
+            double_expanded_nospam_circuits_plus_scratch[i] = double_expanded_ckt
         self.tree = _EvalTree.create(double_expanded_nospam_circuits_plus_scratch)
         #print("Atom tree: %d circuits => tree of size %d" % (len(expanded_nospam_circuits), len(self.tree)))
 
@@ -138,7 +148,7 @@ def add_expanded_circuits(indices, add_to_this_dict):
         # quantity plus a spam-tuple. We order the final indices so that all the outcomes corresponding to a
         # given spam-tuple are contiguous.
 
-        tree_indices_by_spamtuple = _collections.OrderedDict()  # "tree" indices index expanded_nospam_circuits
+        tree_indices_by_spamtuple = dict()  # "tree" indices index expanded_nospam_circuits
         for i, c in expanded_nospam_circuits.items():
             for spam_tuple in expanded_nospam_circuit_outcomes[c].keys():
                 if spam_tuple not in tree_indices_by_spamtuple: tree_indices_by_spamtuple[spam_tuple] = []
@@ -147,7 +157,7 @@ def add_expanded_circuits(indices, add_to_this_dict):
         #Assign element indices, starting at `offset`
         # now that we know how many of each spamtuple there are, assign final element indices.
         local_offset = 0
-        self.indices_by_spamtuple = _collections.OrderedDict()  # values are (element_indices, tree_indices) tuples.
+        self.indices_by_spamtuple = dict()  # values are (element_indices, tree_indices) tuples.
         for spam_tuple, tree_indices in tree_indices_by_spamtuple.items():
             self.indices_by_spamtuple[spam_tuple] = (slice(local_offset, local_offset + len(tree_indices)),
                                                      _slct.list_to_slice(tree_indices, array_ok=True))
@@ -157,8 +167,7 @@ def add_expanded_circuits(indices, add_to_this_dict):
         element_slice = None  # slice(offset, offset + local_offset)  # *global* (of parent layout) element-index slice
         num_elements = local_offset
 
-        elindex_outcome_tuples = _collections.OrderedDict([
-            (unique_i, list()) for unique_i in range(len(unique_complete_circuits))])
+        elindex_outcome_tuples = {unique_i: list() for unique_i in range(len(unique_complete_circuits))}
 
         for spam_tuple, (element_indices, tree_indices) in self.indices_by_spamtuple.items():
             for elindex, tree_index in zip(_slct.indices(element_indices), _slct.to_array(tree_indices)):
@@ -231,7 +240,7 @@ class MatrixCOPALayout(_DistributableCOPALayout):
     Parameters
     ----------
     circuits : list
-        A list of:class:`Circuit` objects representing the circuits this layout will include.
+        A list of :class:`Circuit` objects representing the circuits this layout will include.
 
     model : Model
         The model that will be used to compute circuit outcome probabilities using this layout.
@@ -271,11 +280,16 @@ class MatrixCOPALayout(_DistributableCOPALayout):
     verbosity : int or VerbosityPrinter
         Determines how much output to send to stdout.  0 means no output, higher
         integers mean more output.
+
+    layout_creation_circuit_cache : dict, optional (default None)
+        A precomputed dictionary serving as a cache for completed
+        circuits. I.e. circuits with prep labels and POVM labels appended.
     """
 
     def __init__(self, circuits, model, dataset=None, num_sub_trees=None, num_tree_processors=1,
                  num_param_dimension_processors=(), param_dimensions=(),
-                 param_dimension_blk_sizes=(), resource_alloc=None, verbosity=0):
+                 param_dimension_blk_sizes=(), resource_alloc=None, verbosity=0, 
+                 layout_creation_circuit_cache = None):
 
         #OUTDATED: TODO - revise this:
         # 1. pre-process => get complete circuits => spam-tuples list for each no-spam circuit (no expanding yet)
@@ -290,16 +304,56 @@ def __init__(self, circuits, model, dataset=None, num_sub_trees=None, num_tree_p
         unique_circuits, to_unique = self._compute_unique_circuits(circuits)
         aliases = circuits.op_label_aliases if isinstance(circuits, _CircuitList) else None
         ds_circuits = _lt.apply_aliases_to_circuits(unique_circuits, aliases)
-        unique_complete_circuits, split_unique_circuits = model.complete_circuits(unique_circuits, return_split=True)
+
+        #extract subcaches from layout_creation_circuit_cache:
+        if layout_creation_circuit_cache is None:
+            layout_creation_circuit_cache = dict()
+        self.completed_circuit_cache = layout_creation_circuit_cache.get('completed_circuits', None)
+        self.split_circuit_cache = layout_creation_circuit_cache.get('split_circuits', None)
+        self.expanded_and_separated_circuits_cache = layout_creation_circuit_cache.get('expanded_and_separated_circuits', None)
+        self.expanded_subcircuits_no_spam_cache = layout_creation_circuit_cache.get('expanded_subcircuits_no_spam', None)
+        
+        if self.completed_circuit_cache is None:
+            unique_complete_circuits, split_unique_circuits = model.complete_circuits(unique_circuits, return_split=True)
+        else:
+            unique_complete_circuits = []
+            for c in unique_circuits:
+                comp_ckt = self.completed_circuit_cache.get(c, None)
+                if comp_ckt is not None:
+                    unique_complete_circuits.append(comp_ckt)
+                else:
+                    unique_complete_circuits.append(model.complete_circuit(c))
+
         #Note: "unique" means a unique circuit *before* circuit-completion, so there could be duplicate
         # "unique circuits" after completion, e.g. "rho0Gx" and "Gx" could both complete to "rho0GxMdefault_0".
 
-        circuits_by_unique_nospam_circuits = _collections.OrderedDict()
-        for i, (_, nospam_c, _) in enumerate(split_unique_circuits):
-            if nospam_c in circuits_by_unique_nospam_circuits:
-                circuits_by_unique_nospam_circuits[nospam_c].append(i)
-            else:
-                circuits_by_unique_nospam_circuits[nospam_c] = [i]
+        circuits_by_unique_nospam_circuits = dict()
+        if self.completed_circuit_cache is None:
+            for i, (_, nospam_c, _) in enumerate(split_unique_circuits):
+                if nospam_c in circuits_by_unique_nospam_circuits:
+                    circuits_by_unique_nospam_circuits[nospam_c].append(i)
+                else:
+                    circuits_by_unique_nospam_circuits[nospam_c] = [i]
+            #also create the split circuit cache at this point for future use.
+            if self.split_circuit_cache is None:
+                self.split_circuit_cache = {unique_ckt:split_ckt for unique_ckt, split_ckt in zip(unique_circuits, split_unique_circuits)}
+
+        else:
+            if self.split_circuit_cache is None:
+                self.split_circuit_cache = dict()
+            for i, (c_unique_complete, c_unique) in enumerate(zip(unique_complete_circuits, unique_circuits)):
+                split_ckt_tup = self.split_circuit_cache.get(c_unique, None)
+                nospam_c= split_ckt_tup[1] if split_ckt_tup is not None else None
+                if nospam_c is None:
+                    split_ckt_tup = model.split_circuit(c_unique_complete)
+                    nospam_c= split_ckt_tup[1]
+                    #also add this missing circuit to the cache for future use.
+                    self.split_circuit_cache[c_unique] = split_ckt_tup
+                if nospam_c in circuits_by_unique_nospam_circuits:
+                    circuits_by_unique_nospam_circuits[nospam_c].append(i)
+                else:
+                    circuits_by_unique_nospam_circuits[nospam_c] = [i]
+
         unique_nospam_circuits = list(circuits_by_unique_nospam_circuits.keys())
         
         # Split circuits into groups that will make good subtrees (all procs do this)
@@ -317,9 +371,16 @@ def _create_atom(args):
             group, helpful_scratch_group = args
             return _MatrixCOPALayoutAtom(unique_complete_circuits, unique_nospam_circuits,
                                          circuits_by_unique_nospam_circuits, ds_circuits,
-                                         group, helpful_scratch_group, model, dataset)
+                                         group, helpful_scratch_group, model, 
+                                         unique_circuits, dataset,
+                                         self.expanded_and_separated_circuits_cache,
+                                         self.expanded_subcircuits_no_spam_cache)
 
         super().__init__(circuits, unique_circuits, to_unique, unique_complete_circuits,
                          _create_atom, list(zip(groups, helpful_scratch)), num_tree_processors,
                          num_param_dimension_processors, param_dimensions,
                          param_dimension_blk_sizes, resource_alloc, verbosity)
+
+
+
+
diff --git a/pygsti/layouts/termlayout.py b/pygsti/layouts/termlayout.py
index 10b5458b7..95835c79a 100644
--- a/pygsti/layouts/termlayout.py
+++ b/pygsti/layouts/termlayout.py
@@ -57,7 +57,7 @@ def __init__(self, unique_complete_circuits, ds_circuits, group, model, dataset)
         expanded_circuit_outcomes = _collections.OrderedDict()
         for i in group:
             observed_outcomes = None if (dataset is None) else dataset[ds_circuits[i]].outcomes
-            d = unique_complete_circuits[i].expand_instruments_and_separate_povm(model, observed_outcomes)
+            d = model.expand_instruments_and_separate_povm(unique_complete_circuits[i], observed_outcomes)
             expanded_circuit_outcomes_by_unique[i] = d
             expanded_circuit_outcomes.update(d)
 
diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 368308a46..f3a90f09e 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -32,6 +32,7 @@
 from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
 from pygsti.tools import slicetools as _slct
 from pygsti.tools import matrixtools as _mt
+from pygsti.circuits import Circuit as _Circuit, SeparatePOVMCircuit as _SeparatePOVMCircuit
 
 MEMLIMIT_FOR_NONGAUGE_PARAMS = None
 
@@ -1150,10 +1151,45 @@ def circuit_outcomes(self, circuit):
 
         Returns
         -------
-        tuple
+        tuple corresponding to the possible outcomes for circuit.
         """
-        outcomes = circuit.expand_instruments_and_separate_povm(self)  # dict w/keys=sep-povm-circuits, vals=outcomes
+        outcomes = self.expand_instruments_and_separate_povm(circuit)  # dict w/keys=sep-povm-circuits, vals=outcomes
         return tuple(_itertools.chain(*outcomes.values()))  # concatenate outputs from all sep-povm-circuits
+    
+    def bulk_circuit_outcomes(self, circuits, split_circuits=None, completed_circuits=None):
+        """
+        Get all the possible outcome labels produced by simulating each of the circuits
+        in this list of circuits.
+
+        Parameters
+        ----------
+        circuits : list of Circuits
+            list of Circuits to get outcomes of.
+        
+        split_circuits : list of tuples, optional (default None)
+            If specified, this is a list of tuples for each circuit corresponding to the splitting of
+            the circuit into the prep label, spam-free circuit, and povm label. This is the same format
+            produced by the :meth:split_circuit(s) method, and so this option can allow for accelerating this
+            method when that has previously been run. When using this kwarg only one of this or 
+            the `complete_circuits` kwargs should be used.
+
+        completed_circuits : list of Circuits, optional (default None)
+            If specified, this is a list of compeleted circuits with prep and povm labels included.
+            This is the format produced by the :meth:complete_circuit(s) method, and this can
+            be used to accelerate this method call when that has been previously run. Should not
+            be used in conjunction with `split_circuits`.
+
+        Returns
+        -------
+        list of tuples corresponding to the possible outcomes for each circuit.
+        """
+
+        # list of dict w/keys=sep-povm-circuits, vals=outcomes
+        outcomes_list = self.bulk_expand_instruments_and_separate_povm(circuits, 
+                                                                       split_circuits=split_circuits,
+                                                                       completed_circuits=completed_circuits)  
+        
+        return [tuple(_itertools.chain(*outcomes.values())) for outcomes in outcomes_list]  # concatenate outputs from all sep-povm-circuits
 
     def split_circuit(self, circuit, erroron=('prep', 'povm'), split_prep=True, split_povm=True):
         """
@@ -1356,9 +1392,185 @@ def complete_circuit(self, circuit, prep_lbl_to_prepend=None, povm_lbl_to_append
         Circuit
             Possibly the same object as `circuit`, if no additions are needed.
         """
-
         comp_circuit = self.complete_circuits([circuit], prep_lbl_to_prepend, povm_lbl_to_append, False)
         return comp_circuit[0]
+    
+    def expand_instruments_and_separate_povm(self, circuit, observed_outcomes=None):
+        """
+        Creates a dictionary of :class:`SeparatePOVMCircuit` objects from expanding the instruments of this circuit.
+
+        Each key of the returned dictionary replaces the instruments in this circuit with a selection
+        of their members.  (The size of the resulting dictionary is the product of the sizes of
+        each instrument appearing in this circuit when `observed_outcomes is None`).  Keys are stored
+        as :class:`SeparatePOVMCircuit` objects so it's easy to keep track of which POVM outcomes (effects)
+        correspond to observed data.  This function is, for the most part, used internally to process
+        a circuit before computing its outcome probabilities.
+
+        Parameters
+        ----------
+        circuit : Circuit
+            The circuit to expand, using necessary details regarding the expansion from this model, including:
+
+            - default SPAM layers
+            - definitions of instrument-containing layers
+            - expansions of individual instruments and POVMs
+
+        observed_outcomes : iterable, optional (default None)
+            If specified an iterable over the subset of outcomes empirically observed for this circuit.
+
+        Returns
+        -------
+        OrderedDict
+            A dict whose keys are :class:`SeparatePOVMCircuit` objects and whose
+            values are tuples of the outcome labels corresponding to this circuit,
+            one per POVM effect held in the key.
+        """
+        expanded_circuit_outcomes = self.bulk_expand_instruments_and_separate_povm([circuit], [observed_outcomes])
+        return expanded_circuit_outcomes[0]
+    
+    def bulk_expand_instruments_and_separate_povm(self, circuits, observed_outcomes_list=None, split_circuits = None, 
+                                                  completed_circuits = None):
+        """
+        Creates a list of dictionaries mapping from :class:`SeparatePOVMCircuit` 
+        objects from expanding the instruments of this circuit.
+
+        Each key of the returned dictionary replaces the instruments in this circuit with a selection
+        of their members.  (The size of the resulting dictionary is the product of the sizes of
+        each instrument appearing in this circuit when `observed_outcomes is None`).  Keys are stored
+        as :class:`SeparatePOVMCircuit` objects so it's easy to keep track of which POVM outcomes (effects)
+        correspond to observed data.  This function is, for the most part, used internally to process
+        a circuit before computing its outcome probabilities.
+
+        This function works similarly to expand_instruments_and_separate_povm, except it operates on
+        an entire list of circuits at once, and provides additional kwargs to accelerate computation.
+
+        Parameters
+        ----------
+        circuit : Circuit
+            The circuit to expand, using necessary details regarding the expansion from this model, including:
+
+            - default SPAM layers
+            - definitions of instrument-containing layers
+            - expansions of individual instruments and POVMs
+
+        observed_outcomes_list : list of iterables, optional (default None)
+            If specified a list of iterables over the subset of outcomes empirically observed for each circuit.
+        
+        split_circuits : list of tuples, optional (default None)
+            If specified, this is a list of tuples for each circuit corresponding to the splitting of
+            the circuit into the prep label, spam-free circuit, and povm label. This is the same format
+            produced by the :meth:split_circuit(s) method, and so this option can allow for accelerating this
+            method when that has previously been run. When using this kwarg only one of this or 
+            the `complete_circuits` kwargs should be used.
+
+        completed_circuits : list of Circuits, optional (default None)
+            If specified, this is a list of compeleted circuits with prep and povm labels included.
+            This is the format produced by the :meth:complete_circuit(s) method, and this can
+            be used to accelerate this method call when that has been previously run. Should not
+            be used in conjunction with `split_circuits`.
+
+        Returns
+        -------
+        list of OrderedDicts
+            A list of dictionaries whose keys are :class:`SeparatePOVMCircuit` objects and whose
+            values are tuples of the outcome labels corresponding to each circuit,
+            one per POVM effect held in the key.
+        """
+
+        assert(not (completed_circuits is not None and split_circuits is not None)), "Inclusion of non-trivial values"\
+              +" for both `complete_circuits` and `split_circuits` is not supported. Please use only one of these two arguments."
+
+        if split_circuits is not None:
+            povm_lbls = [split_ckt[2] for split_ckt in split_circuits]
+            circuits_without_povm = [(split_ckt[0],) + split_ckt[1] for split_ckt in split_circuits]
+        elif completed_circuits is not None:
+            povm_lbls = [comp_ckt[-1] for comp_ckt in completed_circuits]
+            circuits_without_povm = [comp_ckt[:-1] for comp_ckt in completed_circuits]
+        else:
+            completed_circuits = self.complete_circuits(circuits)
+            povm_lbls = [comp_ckt[-1] for comp_ckt in completed_circuits]
+            circuits_without_povm = [comp_ckt[:-1] for comp_ckt in completed_circuits]
+        
+        if observed_outcomes_list is None:
+            observed_outcomes_list = [None]*len(circuits)
+
+
+        expanded_circuit_outcomes_list = [_collections.OrderedDict() for _ in range(len(circuits))]
+
+        def create_tree(lst):
+            subs = _collections.OrderedDict()
+            for el in lst:
+                if len(el) > 0:
+                    if el[0] not in subs: subs[el[0]] = []
+                    subs[el[0]].append(el[1:])
+            return _collections.OrderedDict([(k, create_tree(sub_lst)) for k, sub_lst in subs.items()])
+
+        def add_expanded_circuit_outcomes(circuit, running_outcomes, ootree, start):
+            """
+            """
+            cir = circuit if start == 0 else circuit[start:]  # for performance, avoid uneeded slicing
+            for k, layer_label in enumerate(cir, start=start):
+                components = layer_label.components
+                #instrument_inds = _np.nonzero([model._is_primitive_instrument_layer_lbl(component)
+                #                               for component in components])[0]  # SLOWER than statement below
+                instrument_inds = _np.array([i for i, component in enumerate(components)
+                                             if self._is_primitive_instrument_layer_lbl(component)])
+                if instrument_inds.size > 0:
+                    # This layer contains at least one instrument => recurse with instrument(s) replaced with
+                    #  all combinations of their members.
+                    component_lookup = {i: comp for i, comp in enumerate(components)}
+                    instrument_members = [self._member_labels_for_instrument(components[i])
+                                          for i in instrument_inds]  # also components of outcome labels
+                    for selected_instrmt_members in _itertools.product(*instrument_members):
+                        expanded_layer_lbl = component_lookup.copy()
+                        expanded_layer_lbl.update({i: components[i] + "_" + sel
+                                                   for i, sel in zip(instrument_inds, selected_instrmt_members)})
+                        expanded_layer_lbl = _Label([expanded_layer_lbl[i] for i in range(len(components))])
+
+                        if ootree is not None:
+                            new_ootree = ootree
+                            for sel in selected_instrmt_members:
+                                new_ootree = new_ootree.get(sel, {})
+                            if len(new_ootree) == 0: continue  # no observed outcomes along this outcome-tree path
+                        else:
+                            new_ootree = None
+
+                        add_expanded_circuit_outcomes(circuit[0:k] + _Circuit((expanded_layer_lbl,)) + circuit[k + 1:],
+                                                      running_outcomes + selected_instrmt_members, new_ootree, k + 1)
+                    break
+
+            else:  # no more instruments to process: `cir` contains no instruments => add an expanded circuit
+                assert(circuit not in expanded_circuit_outcomes)  # shouldn't be possible to generate duplicates...
+                elabels = self._effect_labels_for_povm(povm_lbl) if (observed_outcomes is None) \
+                    else tuple(ootree.keys())
+                outcomes = tuple((running_outcomes + (elabel,) for elabel in elabels))
+                expanded_circuit_outcomes[_SeparatePOVMCircuit(circuit, povm_lbl, elabels)] = outcomes
+
+        has_instruments = self._has_instruments()
+        unique_povm_labels = set(povm_lbls)
+        effect_label_dict = {povm_lbl: self._effect_labels_for_povm(povm_lbl) for povm_lbl in unique_povm_labels}
+
+        for povm_lbl, circuit_without_povm, expanded_circuit_outcomes, observed_outcomes in zip(povm_lbls, circuits_without_povm, 
+                                                                                                expanded_circuit_outcomes_list, 
+                                                                                                observed_outcomes_list):
+            ootree = create_tree(observed_outcomes) if observed_outcomes is not None else None  # tree of observed outcomes
+            # e.g. [('0','00'), ('0','01'), ('1','10')] ==> {'0': {'00': {}, '01': {}}, '1': {'10': {}}}
+
+            if has_instruments:
+                add_expanded_circuit_outcomes(circuit_without_povm, (), ootree, start=0)
+            else:
+                # It may be helpful to cache the set of elabels for a POVM (maybe within the model?) because
+                # currently the call to _effect_labels_for_povm may be a bottleneck.  It's needed, even when we have
+                # observed outcomes, because there may be some observed outcomes that aren't modeled (e.g. leakage states)
+                if observed_outcomes is None:
+                    elabels = effect_label_dict[povm_lbl]
+                else:
+                    possible_lbls = set(effect_label_dict[povm_lbl])
+                    elabels = tuple([oo for oo in ootree.keys() if oo in possible_lbls])
+                outcomes = tuple(((elabel,) for elabel in elabels))
+                expanded_circuit_outcomes[_SeparatePOVMCircuit(circuit_without_povm, povm_lbl, elabels)] = outcomes
+
+        return expanded_circuit_outcomes_list
 
     def complete_circuits(self, circuits, prep_lbl_to_prepend=None, povm_lbl_to_append=None, return_split = False):
         """
diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py
index 191fd736b..208bdb46d 100644
--- a/pygsti/objectivefns/objectivefns.py
+++ b/pygsti/objectivefns/objectivefns.py
@@ -19,11 +19,13 @@
 
 from pygsti import tools as _tools
 from pygsti.layouts.distlayout import DistributableCOPALayout as _DistributableCOPALayout
+from pygsti.layouts.matrixlayout import MatrixCOPALayout as _MatrixCOPALayout
 from pygsti.tools import slicetools as _slct, mpitools as _mpit, sharedmemtools as _smt
 from pygsti.circuits.circuitlist import CircuitList as _CircuitList
 from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation
 from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 from pygsti.baseobjs.verbosityprinter import VerbosityPrinter as _VerbosityPrinter
+from pygsti.models.model import OpModel as _OpModel
 
 
 def _objfn(objfn_cls, model, dataset, circuits=None,
@@ -843,12 +845,10 @@ class ModelDatasetCircuitsStore(object):
     point.
     """
     def __init__(self, model, dataset, circuits=None, resource_alloc=None, array_types=(),
-                 precomp_layout=None, verbosity=0):
+                 precomp_layout=None, verbosity=0, outcome_count_by_circuit=None):
         self.dataset = dataset
         self.model = model
-        #self.opBasis = mdl.basis
         self.resource_alloc = _ResourceAllocation.cast(resource_alloc)
-        # expand = ??? get from model based on fwdsim type?
 
         circuit_list = circuits if (circuits is not None) else list(dataset.keys())
         bulk_circuit_list = circuit_list if isinstance(
@@ -872,8 +872,21 @@ def __init__(self, model, dataset, circuits=None, resource_alloc=None, array_typ
         else:
             self.global_circuits = self.circuits
 
-        #self.circuits = bulk_circuit_list[:]
-        #self.circuit_weights = bulk_circuit_list.circuit_weights
+        #If a matrix layout then we have some precached circuit structures we can
+        #grab to speed up store generation.
+        if isinstance(self.layout, _MatrixCOPALayout):
+            #Grab the split_circuit_cache and down select to those in
+            #self.circuits
+            self.split_circuit_cache = self.layout.split_circuit_cache
+            self.split_circuits = [self.split_circuit_cache[ckt] for ckt in self.circuits]
+        #currently only implemented for matrix, will eventually add map support.
+        else:
+            self.split_circuits = None
+            self.split_circuit_cache = None
+
+        #set the value of the circuit outcome count cache (can be None)
+        self.outcome_count_by_circuit_cache = outcome_count_by_circuit
+
         self.ds_circuits = self.circuits.apply_aliases()
 
         # computed by add_count_vectors
@@ -888,18 +901,6 @@ def __init__(self, model, dataset, circuits=None, resource_alloc=None, array_typ
 
         self.time_dependent = False  # indicates whether the data should be treated as time-resolved
 
-        #if not self.cache.has_evaltree():
-        #    subcalls = self.get_evaltree_subcalls()
-        #    evt_resource_alloc = _ResourceAllocation(self.raw_objfn.comm, evt_mlim,
-        #                                             self.raw_objfn.profiler, self.raw_objfn.distribute_method)
-        #    self.cache.add_evaltree(self.mdl, self.dataset, bulk_circuit_list, evt_resource_alloc,
-        #                            subcalls, self.raw_objfn.printer - 1)
-        #self.eval_tree = self.cache.eval_tree
-        #self.lookup = self.cache.lookup
-        #self.outcomes_lookup = self.cache.outcomes_lookup
-        #self.wrt_block_size = self.cache.wrt_block_size
-        #self.wrt_block_size2 = self.cache.wrt_block_size2
-
         #convenience attributes (could make properties?)
         if isinstance(self.layout, _DistributableCOPALayout):
             self.global_nelements = self.layout.global_num_elements
@@ -940,14 +941,36 @@ def add_omitted_freqs(self, printer=None, force=False):
         """
         if self.firsts is None or force:
             # FUTURE: add any tracked memory? self.resource_alloc.add_tracked_memory(...)
-            self.firsts = []; self.indicesOfCircuitsWithOmittedData = []
-            for i, c in enumerate(self.circuits):
-                indices = _slct.to_array(self.layout.indices_for_index(i))
-                lklen = _slct.length(self.layout.indices_for_index(i))
-                if 0 < lklen < self.model.compute_num_outcomes(c):
+            self.firsts = [] 
+            self.indicesOfCircuitsWithOmittedData = []
+
+            if self.outcome_count_by_circuit_cache is None:
+                #bulk compute the number of outcomes.
+                if isinstance(self.model, _OpModel) and self.split_circuits is not None:
+                    bulk_outcomes_list = self.model.bulk_circuit_outcomes(self.circuits, split_circuits=self.split_circuits)
+                    num_outcomes_list = [len(outcome_tup) for outcome_tup in bulk_outcomes_list]
+                else:
+                    num_outcomes_list = [self.model.compute_num_outcomes(c) for c in self.circuits]
+            else:
+                num_outcomes_list = []
+                for ckt in self.circuits:
+                    num_outcomes = self.outcome_count_by_circuit_cache.get(ckt, None)
+                    if num_outcomes is None:
+                        num_outcomes = self.model.compute_num_outcomes(ckt)
+                        #also add this to the cache, just in case it is later needed.
+                        self.outcome_count_by_circuit_cache[ckt] = num_outcomes
+                    num_outcomes_list.append(num_outcomes)
+
+            for i in range(len(self.circuits)):
+                indices = self.layout.indices_for_index(i)
+                #The return types of indices_for_index are either ndarrays
+                #or slices.
+                if isinstance(indices, slice): 
+                    indices = _slct.indices(indices)
+                if 0 < len(indices) < num_outcomes_list[i]:
                     self.firsts.append(indices[0])
                     self.indicesOfCircuitsWithOmittedData.append(i)
-            if len(self.firsts) > 0:
+            if self.firsts:
                 self.firsts = _np.array(self.firsts, 'i')
                 self.indicesOfCircuitsWithOmittedData = _np.array(self.indicesOfCircuitsWithOmittedData, 'i')
                 self.dprobs_omitted_rowsum = _np.empty((len(self.firsts), self.nparams), 'd')
@@ -974,13 +997,15 @@ def add_count_vectors(self, force=False):
 
             for (i, circuit) in enumerate(self.ds_circuits):
                 cnts = self.dataset[circuit].counts
-                totals[self.layout.indices_for_index(i)] = sum(cnts.values())  # dataset[opStr].total
-                counts[self.layout.indices_for_index(i)] = [cnts.get(x, 0) for x in self.layout.outcomes_for_index(i)]
+                idcs_for_idx = self.layout.indices_for_index(i)
+                totals[idcs_for_idx] = sum(cnts.values())  # dataset[opStr].
+                counts[idcs_for_idx] = [cnts.getitem_unsafe(x, 0) for x in self.layout.outcomes_for_index(i)]
 
             if self.circuits.circuit_weights is not None:
                 for i in range(len(self.ds_circuits)):  # multiply N's by weights
-                    counts[self.layout.indices_for_index(i)] *= self.circuits.circuit_weights[i]
-                    totals[self.layout.indices_for_index(i)] *= self.circuits.circuit_weights[i]
+                    idcs_for_idx = self.layout.indices_for_index(i)
+                    counts[idcs_for_idx] *= self.circuits.circuit_weights[i]
+                    totals[idcs_for_idx] *= self.circuits.circuit_weights[i]
 
             self.counts = counts
             self.total_counts = totals
@@ -994,7 +1019,7 @@ class EvaluatedModelDatasetCircuitsStore(ModelDatasetCircuitsStore):
 
     def __init__(self, mdc_store, verbosity):
         super().__init__(mdc_store.model, mdc_store.dataset, mdc_store.global_circuits, mdc_store.resource_alloc,
-                         mdc_store.array_types, mdc_store.layout, verbosity)
+                         mdc_store.array_types, mdc_store.layout, verbosity, mdc_store.outcome_count_by_circuit_cache)
 
         # Memory check - see if there's enough memory to hold all the evaluated quantities
         #persistent_mem = self.layout.memory_estimate()
diff --git a/pygsti/tools/slicetools.py b/pygsti/tools/slicetools.py
index ba49b1056..506045182 100644
--- a/pygsti/tools/slicetools.py
+++ b/pygsti/tools/slicetools.py
@@ -26,7 +26,8 @@ def length(s):
     -------
     int
     """
-    if not isinstance(s, slice): return len(s)
+    if not isinstance(s, slice): 
+        return len(s)
     if s.start is None or s.stop is None:
         return 0
     if s.step is None:
@@ -191,7 +192,8 @@ def indices(s, n=None):
     elif s.start < 0:
         assert(n is not None), "Must supply `n` to obtain indices of a slice with negative start point!"
         start = n + s.start
-    else: start = s.start
+    else: 
+        start = s.start
 
     if s.stop is None:
         assert(n is not None), "Must supply `n` to obtain indices of a slice with unspecified stop point!"
@@ -199,12 +201,56 @@ def indices(s, n=None):
     elif s.stop < 0:
         assert(n is not None), "Must supply `n` to obtain indices of a slice with negative stop point!"
         stop = n + s.stop
-    else: stop = s.stop
+    else: 
+        stop = s.stop
 
     if s.step is None:
         return list(range(start, stop))
-    return list(range(start, stop, s.step))
+    else:
+        return list(range(start, stop, s.step))
+    
+def indices_as_array(s, n=None):
+    """
+    Returns a numpy array of the indices specified by slice `s`.
+
+    Parameters
+    ----------
+    s : slice
+        The slice to operate upon.
+
+    n : int, optional
+        The number of elements in the array being indexed,
+        used for computing *negative* start/stop points.
+
+    Returns
+    -------
+    numpy ndarray array of integers
+    """
+    if s.start is None and s.stop is None:
+        return []
+
+    if s.start is None:
+        start = 0
+    elif s.start < 0:
+        assert(n is not None), "Must supply `n` to obtain indices of a slice with negative start point!"
+        start = n + s.start
+    else: 
+        start = s.start
+
+    if s.stop is None:
+        assert(n is not None), "Must supply `n` to obtain indices of a slice with unspecified stop point!"
+        stop = n
+    elif s.stop < 0:
+        assert(n is not None), "Must supply `n` to obtain indices of a slice with negative stop point!"
+        stop = n + s.stop
+    else: 
+        stop = s.stop
 
+    if s.step is None:
+        return _np.arange(start, stop, dtype=_np.int64)
+    else:
+        return _np.arange(start, stop, s.step, dtype=_np.int64)
+    
 
 def list_to_slice(lst, array_ok=False, require_contiguous=True):
     """
@@ -240,17 +286,23 @@ def list_to_slice(lst, array_ok=False, require_contiguous=True):
                 else:
                     raise ValueError("Slice must be contiguous!")
         return lst
-    if lst is None or len(lst) == 0: return slice(0, 0)
+    if lst is None or len(lst) == 0: 
+        return slice(0, 0)
     start = lst[0]
 
-    if len(lst) == 1: return slice(start, start + 1)
-    step = lst[1] - lst[0]; stop = start + step * len(lst)
+    if len(lst) == 1: 
+        return slice(start, start + 1)
+    step = lst[1] - lst[0]
+    stop = start + step * len(lst)
 
     if list(lst) == list(range(start, stop, step)):
         if require_contiguous and step != 1:
-            if array_ok: return _np.array(lst, _np.int64)
-            else: raise ValueError("Slice must be contiguous (or array_ok must be True)!")
-        if step == 1: step = None
+            if array_ok: 
+                return _np.array(lst, _np.int64)
+            else: 
+                raise ValueError("Slice must be contiguous (or array_ok must be True)!")
+        if step == 1: 
+            step = None
         return slice(start, stop, step)
     elif array_ok:
         return _np.array(lst, _np.int64)
@@ -272,7 +324,7 @@ def to_array(slc_or_list_like):
     numpy.ndarray
     """
     if isinstance(slc_or_list_like, slice):
-        return _np.array(indices(slc_or_list_like), _np.int64)
+        return indices_as_array(slc_or_list_like)
     else:
         return _np.array(slc_or_list_like, _np.int64)
 


Test patch:
diff --git a/test/unit/protocols/test_rb.py b/test/unit/protocols/test_rb.py
index fcb616cf2..9f461f098 100644
--- a/test/unit/protocols/test_rb.py
+++ b/test/unit/protocols/test_rb.py
@@ -1,6 +1,7 @@
 from ..util import BaseCase
 
 import numpy as _np
+from pathlib import Path
 
 import pygsti
 from pygsti.protocols import rb as _rb
@@ -9,6 +10,8 @@
 from pygsti.circuits import Circuit
 from pygsti.baseobjs import Label
 
+FILE_PATH = str(Path(__file__).resolve().parent)
+
 class TestCliffordRBDesign(BaseCase):
 
     def setUp(self):
@@ -111,9 +114,9 @@ def test_serialization(self):
             citerations=self.citerations, compilerargs=self.compiler_args, seed=self.seed,
             verbosity=self.verbosity, num_processes=1)
 
-        crb_design.write('../../test_packages/temp_test_files/test_CliffordRBDesign_serialization')
+        crb_design.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_CliffordRBDesign_serialization')
         #then read this back in
-        crb_design_read = _rb.CliffordRBDesign.from_dir('../../test_packages/temp_test_files/test_CliffordRBDesign_serialization')
+        crb_design_read = _rb.CliffordRBDesign.from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_CliffordRBDesign_serialization')
 
         self.assertEqual(crb_design.all_circuits_needing_data, crb_design_read.all_circuits_needing_data)
         self.assertEqual(crb_design.interleaved_circuit, crb_design_read.interleaved_circuit)
@@ -163,9 +166,9 @@ def test_combined_design_access(self):
 
     def test_serialization(self):
 
-        self.irb_design.write('../../test_packages/temp_test_files/test_InterleavedRBDesign_serialization')
+        self.irb_design.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_InterleavedRBDesign_serialization')
         #then read this back in
-        irb_design_read = _rb.InterleavedRBDesign.from_dir('../../test_packages/temp_test_files/test_InterleavedRBDesign_serialization')
+        irb_design_read = _rb.InterleavedRBDesign.from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_InterleavedRBDesign_serialization')
 
         self.assertEqual(self.irb_design.all_circuits_needing_data, irb_design_read.all_circuits_needing_data)
         self.assertEqual(self.irb_design['crb'].all_circuits_needing_data, irb_design_read['crb'].all_circuits_needing_data)
@@ -248,9 +251,9 @@ def test_serialization(self):
             conditionaltwirl=True, citerations=self.citerations, compilerargs=self.compiler_args,
             partitioned=False, seed=self.seed, verbosity=self.verbosity, num_processes=1)
 
-        drb_design.write('../../test_packages/temp_test_files/test_DirectRBDesign_serialization')
+        drb_design.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_DirectRBDesign_serialization')
         #then read this back in
-        drb_design_read = _rb.DirectRBDesign.from_dir('../../test_packages/temp_test_files/test_DirectRBDesign_serialization')
+        drb_design_read = _rb.DirectRBDesign.from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_DirectRBDesign_serialization')
 
         self.assertEqual(drb_design.all_circuits_needing_data, drb_design_read.all_circuits_needing_data)
 
@@ -375,9 +378,9 @@ def test_serialization(self):
             localclifford=True, paulirandomize=True, seed=self.seed, verbosity=self.verbosity,
             num_processes=1)
 
-        mrb_design.write('../../test_packages/temp_test_files/test_MirrorRBDesign_serialization')
+        mrb_design.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_MirrorRBDesign_serialization')
         #then read this back in
-        mrb_design_read = _rb.MirrorRBDesign.from_dir('../../test_packages/temp_test_files/test_MirrorRBDesign_serialization')
+        mrb_design_read = _rb.MirrorRBDesign.from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_MirrorRBDesign_serialization')
 
         self.assertEqual(mrb_design.all_circuits_needing_data, mrb_design_read.all_circuits_needing_data)
 
@@ -424,9 +427,9 @@ def test_serialization(self):
                                                  sampler=self.sampler, samplerargs=self.samplerargs, 
                                                  seed=self.seed, verbosity=0)
         
-        birb_design.write('../../test_packages/temp_test_files/test_BinaryRBDesign_serialization')
+        birb_design.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_BinaryRBDesign_serialization')
         #then read this back in
-        birb_design_read = _rb.BinaryRBDesign.from_dir('../../test_packages/temp_test_files/test_BinaryRBDesign_serialization')
+        birb_design_read = _rb.BinaryRBDesign.from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_BinaryRBDesign_serialization')
 
         self.assertEqual(birb_design.all_circuits_needing_data, birb_design_read.all_circuits_needing_data)
         
@@ -533,8 +536,8 @@ def test_cliffordrb_protocol_ideal(self):
         self.assertTrue(abs(result.fits['A-fixed'].estimates['r'])<=3e-5)
 
         #also test writing and reading the results from disk.
-        result.write('../../test_packages/temp_test_files/test_RandomizedBenchmarking_results')
-        result_read = pygsti.io.read_results_from_dir('../../test_packages/temp_test_files/test_RandomizedBenchmarking_results')
+        result.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_RandomizedBenchmarking_results')
+        result_read = pygsti.io.read_results_from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_RandomizedBenchmarking_results')
         
     def test_cliffordrb_protocol_noisy(self):
         proto = pygsti.protocols.rb.RandomizedBenchmarking(datatype='success_probabilities', defaultfit='A-fixed', rtype='EI',
@@ -703,8 +706,8 @@ def test_interleavedrb_protocol_ideal(self):
         self.assertTrue(abs(estimated_irb_num) <= 1e-5)
 
         #also test writing and reading the results from disk.
-        result.write('../../test_packages/temp_test_files/test_InterleavedRandomizedBenchmarking_results')
-        result_read = pygsti.io.read_results_from_dir('../../test_packages/temp_test_files/test_InterleavedRandomizedBenchmarking_results')
+        result.write(f'{FILE_PATH}/../../test_packages/temp_test_files/test_InterleavedRandomizedBenchmarking_results')
+        result_read = pygsti.io.read_results_from_dir(f'{FILE_PATH}/../../test_packages/temp_test_files/test_InterleavedRandomizedBenchmarking_results')
         
         
     def test_interleavedrb_protocol_noisy(self):