Initial draft of instrument support.

[cgranade]

As discussed in #1673, this PR is an initial draft of what instrument support could look like to get feedback and discussion going on the API and functionality. As such, code style, documentation, unit testing and so forth still need to be addressed (see check list below).

Checklist
Thank you for contributing to QuTiP! Please make sure you have finished the following tasks before opening the PR.

• Please read Contributing to QuTiP Development
• Contributions to qutip should follow the pep8 style.
  You can use pycodestyle to check your code automatically
• Please add tests to cover your changes if applicable.
• If the behavior of the code has changed or new feature has been added, please also update the documentation in the doc folder, and the notebook. Feel free to ask if you are not sure.

Delete this checklist after you have completed all the tasks. If you have not finished them all, you can also open a Draft Pull Request to let the others know this on-going work and keep this checklist in the PR description.

Description
This PR adds a new class, QInstrument, that wraps a decomposition of a quantum instrument into completely positive trace non-increasing maps (subnormalized channels). This new class can be used to predict measurement outcomes and post-measurement states for a variety of different quantum systems:

>>> import qutip
>>> H = qutip.operations.hadamard_transform()
>>> ket_plus = H * qutip.basis(2, 0)
>>> z_instrument = qutip.QInstrument.basis_measurement(2)
>>> (H * z_instrument)(ket_plus)
{Seq(0,): Outcome(probability=0.5000000000000002, output_state=Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[0.5 0.5]
 [0.5 0.5]]), Seq(1,): Outcome(probability=0.5000000000000002, output_state=Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[ 0.5 -0.5]
 [-0.5  0.5]])}
>>> z_instrument.sample(ket_plus)
(Seq(1,), Outcome(probability=0.5000000000000001, output_state=Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True
Qobj data =
[[0. 0.]
 [0. 1.]]))

Instruments can be combined through composition and tensor products (corresponding to Seq and Par outcome labels, respectively):

>>> z_instrument * z_instrument
QInstrument id=2e7d878ab80 {
    dims [[[2], [2]], [[2], [2]]]
    outcomes Seq(0, 0) Seq(1, 1)
}
>>> qutip.tensor([z_instrument] * 2)
QInstrument id=2e7d8a1dcd0 {
    dims [[[2, 2], [2, 2]], [[2, 2], [2, 2]]]
    outcomes Par(0, 0) Par(0, 1) Par(1, 0) Par(1, 1)
}

Impossible outcomes are detected and truncated automatically:

>>> z_instrument * z_instrument
QInstrument id=2e7d8a6e130 {
    dims [[[2], [2]], [[2], [2]]]
    outcomes Seq(0, 0) Seq(1, 1)
}
>>> z_instrument.with_finite_visibility(0.95) ** 2
QInstrument id=2e7d8a6ef10 {
    dims [[[2], [2]], [[2], [2]]]
    outcomes Seq(0, 0) Seq(1, 0) Seq(0, 1) Seq(1, 1)
}

Arbitrary subsystem dims are supported:

>>> qutip.QInstrument.basis_measurement(3)
QInstrument id=2e7d5ca4eb0 {
    dims [[[3], [3]], [[3], [3]]]
    outcomes Seq(0,) Seq(1,) Seq(2,)
}

Incomplete instruments (that is, where the probability of obtaining any result is less than 1, as in the erasure channel case) can be completed:

>>> qutip.QInstrument(qutip.projection(4, 0, 0)).complete()
QInstrument id=1f9735f3f10 {
    dims [[[4], [4]], [[4], [4]]]
    outcomes Seq() Seq('⊥',)
}

Measurement outcome labels can be re-indexed onto integer labels according to lexographical sorting of original labels:

>>> qutip.QInstrument(qutip.projection(4, 0, 0)).complete().reindex()
QInstrument id=1f96e774970 {
    dims [[[4], [4]], [[4], [4]]]
    outcomes Seq(0,) Seq(1,)
}

Composition and tensor products can be combined:

>>> qutip.tensor(z_instrument.with_finite_visibility(0.95) ** 2, z_instrument)
QInstrument id=2e7d89a19a0 {
    dims [[[2, 2], [2, 2]], [[2, 2], [2, 2]]]
    outcomes Par(Seq(0, 0), 0) Par(Seq(0, 0), 1) Par(Seq(1, 0), 0) Par(Seq(1, 0), 1) Par(Seq(0, 1), 0) Par(Seq(0, 1), 1) Par(Seq(1, 1), 0) Par(Seq(1, 1), 1)
}

Outcomes can be labeled by arbitrary hashable types:

>>> qutip.QInstrument.pauli_measurement("ZZ") * qutip.QInstrument.pauli_measurement("XX")
QInstrument id=2e7d89d5dc0 {
    dims [[[2, 2], [2, 2]], [[2, 2], [2, 2]]]
    outcomes Seq(+XX, +ZZ) Seq(-XX, +ZZ) Seq(+XX, -ZZ) Seq(-XX, -ZZ)
}
>>> (qutip.QInstrument.pauli_measurement("ZZ") * qutip.QInstrument.pauli_measurement("XX")) ** 2
QInstrument id=2e7d89a44c0 {
    dims [[[2, 2], [2, 2]], [[2, 2], [2, 2]]]
    outcomes Seq(+XX, +ZZ, +XX, +ZZ) Seq(-XX, +ZZ, -XX, +ZZ) Seq(+XX, -ZZ, +XX, -ZZ) Seq(-XX, -ZZ, -XX, -ZZ)
}
>>> (qutip.QInstrument.pauli_measurement("ZZ") * qutip.QInstrument.pauli_measurement("XX")).with_finite_visibility(0.95) ** 2
QInstrument id=2e7d6654520 {
    dims [[[2, 2], [2, 2]], [[2, 2], [2, 2]]]
    outcomes Seq(+XX, +ZZ, +XX, +ZZ) Seq(-XX, +ZZ, +XX, +ZZ) Seq(+XX, -ZZ, +XX, +ZZ) Seq(-XX, -ZZ, +XX, +ZZ) Seq(+XX, +ZZ, -XX, +ZZ) Seq(-XX, +ZZ, -XX, +ZZ) Seq(+XX, -ZZ, -XX, +ZZ) Seq(-XX, -ZZ, 
-XX, +ZZ) Seq(+XX, +ZZ, +XX, -ZZ) Seq(-XX, +ZZ, +XX, -ZZ) Seq(+XX, -ZZ, +XX, -ZZ) Seq(-XX, -ZZ, +XX, -ZZ) Seq(+XX, +ZZ, -XX, -ZZ) Seq(-XX, +ZZ, -XX, -ZZ) Seq(+XX, -ZZ, -XX, -ZZ) Seq(-XX, -ZZ, -XX, -ZZ)
}

Nonselective channels (that is, the channel resulting from discarding all measurement outcomes) can be efficiently obtained using the nonselective_process property:

>>> ins = (qutip.QInstrument.pauli_measurement("ZZ") * qutip.QInstrument.pauli_measurement("XX")).with_finite_visibility(0.95) ** 2
>>> ins.nonselective_process
Quantum object: dims = [[[2, 2], [2, 2]], [[2, 2], [2, 2]]], shape = (16, 16), type = super, isherm = True
Qobj data =
[[0.50125 0.      0.      0.      0.      0.      0.      0.      0.
  0.      0.      0.      0.      0.      0.      0.49875]
 [0.      0.0025  0.      0.      0.      0.      0.      0.      0.
  0.      0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.0025  0.      0.      0.      0.      0.      0.
  0.      0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.50125 0.      0.      0.      0.      0.
  0.      0.      0.      0.49875 0.      0.      0.     ]
 [0.      0.      0.      0.      0.0025  0.      0.      0.      0.
  0.      0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.50125 0.      0.      0.
  0.      0.49875 0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.50125 0.      0.
  0.49875 0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.      0.0025  0.
  0.      0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.      0.      0.0025
  0.      0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.49875 0.      0.
  0.50125 0.      0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.49875 0.      0.      0.
  0.      0.50125 0.      0.      0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.      0.      0.
  0.      0.      0.0025  0.      0.      0.      0.     ]
 [0.      0.      0.      0.49875 0.      0.      0.      0.      0.
  0.      0.      0.      0.50125 0.      0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.      0.      0.
  0.      0.      0.      0.      0.0025  0.      0.     ]
 [0.      0.      0.      0.      0.      0.      0.      0.      0.
  0.      0.      0.      0.      0.      0.0025  0.     ]
 [0.49875 0.      0.      0.      0.      0.      0.      0.      0.
  0.      0.      0.      0.      0.      0.      0.50125]]

Trace preserving and complete positivity conditions can be checked in a similar way to Qobj instances:

>>> ins.iscptp
True
>>> ins.ishp
True
>>> ins.istp
True
>>> ins.iscp
True

Finally, for convienence, this PR also proposes QuaEC-style shorthand for tensor products of quantum objects and instruments:

>>> z_instrument & z_instrument.with_finite_visibility(0.95)
QInstrument id=2e7d5ca4d30 {
    dims [[[2, 2], [2, 2]], [[2, 2], [2, 2]]]
    outcomes Par(0, 0) Par(0, 1) Par(1, 0) Par(1, 1)
}
>>> z_instrument ^ 3
QInstrument id=2e7d8017e80 {
    dims [[[2, 2, 2], [2, 2, 2]], [[2, 2, 2], [2, 2, 2]]]
    outcomes Par(0, 0, 0) Par(0, 0, 1) Par(0, 1, 0) Par(0, 1, 1) Par(1, 0, 0) Par(1, 0, 1) Par(1, 1, 0) Par(1, 1, 1)
}

Related issues or PRs

• Instrument and operator-sum support #1673

Changelog
Adds support for quantum instruments, a generalization of measurements and channels.

[AGaliciaMartinez]

This will catch arbitrary exceptions. It should be more specific to catch only the required exception.

[cgranade]

This is a draft PR to get initial feedback on the API design; it's not intended that code in this PR is in its final working state (e.g. there's not much in the way of unit testing, the use of the walrus operator isn't compatible with some versions of Python that should be supported, and so forth).

That said, the intent of this particular function is to check the iterability of a given value without causing any other effects on execution — it would be surprising for something analogous to an isinstance call to raise an exception, such that exceptions should not propagate out from this check.

[AGaliciaMartinez]

Hey @cgranade, if you planned to address this later, once the discussion on the API was finished, please, feel free to ignore this comment. I was just interested in the PR and decided to take a quick look and, since I saw the catch-all exception, I thought I could point it out to avoid it appearing in the final working version.

The reason I think it is dangerous to catch all in the try/except statement is that it could lead to undesired and inconsistent behaviour, difficult to debug. For instance, a catch-all try/except will also intercept a KeyboardInterrupt which in turn is used to stop the execution of a program (for example in a notebook). The following code when executed in a notebook will not stop when trying to interrupt the kernel:

def _is_iterable(obj):
    try:
        time.sleep(1)
        iter(obj)
        return True
    except:
        return False
while True:
    print(_is_iterable([1]))

Such problem seems unlikely to happen, after all, I required to include a sleep time of 1 second to reproduce the bug. However, it could happen. Such problems can be avoided by catching TypeError instead of any exception.

[hodgestar]

Probably good for us to leave this note here in order to remind us to fix this later, even if it's not important right this second. It would be bad to forget. :)

P.S. We probably only want to catch TypeError or something like that.

[hodgestar]

@cgranade I'm really liking the overall look of the Instrument class -- it's a much easier way to deal with measurements than calling the measurement operations all the time & of combine measurement operations with other operations on Qobjs.

I'm keen to hear feedback from others too, but in the mean time I'm going to note some thoughts here for when I come back to this:

• It would be good to think about how all the operations on Qobj match up with QuTiP version 5 (since this draft is based on v4 currently).
• I'd like to think about removing Seq and Par and replacing them with some simple rules for sequences, strings and numbers. This would match, e.g., qutip.ket("01") and qutip.basis([2, 2], [0, 1]). Seq and Par do however make it really clear that in one case measurements follow each other on the same subspaces and in the other they are performed simultaneously on different subspaces, so I'm not quite sure. Maybe there is some middle ground.
• In QuTiP it's more normal to have the helper constructors not on the class as classmethods. I do like secondary constructors as classmethods, but I think in cases like this where there are essentially an infinite number of possible constructors, it makes sense to not "bless" any of them by sticking them on the class.
• It would be good to have one really nice use case that we could turn into a guide entry in the documentation. The current small examples are great, but it would be good to add one slightly bigger worked example that did something more "exciting".
• I would not push everything into the qutip namespace (largely because we would likely not want that in v5, although I should actually check what v5 is doing in qutip.__init__ these days).
• We should decide whether to target QuTiP 4.7 or 4.7 and 5 for this. Target just v5 means not having to worry about making it nice in both, but will mean it'll be a bit more of a delay before release).

Hoping to hear comments from others!

[hodgestar]

@cgranade When we initially talked about this, we also spoke about storing Kraus superoperators with it, but it's not completely clear to me whether this route covers that case. "No" is fine (we don't have to make this work for everything) but if the answer is "Yes" or "Maybe", what would that look like?

[cgranade]

@cgranade I'm really liking the overall look of the Instrument class -- it's a much easier way to deal with measurements than calling the measurement operations all the time & of combine measurement operations with other operations on Qobjs.

Thank you for the kind words!

I'm keen to hear feedback from others too, but in the mean time I'm going to note some thoughts here for when I come back to this:

• It would be good to think about how all the operations on Qobj match up with QuTiP version 5 (since this draft is based on v4 currently).
• We should decide whether to target QuTiP 4.7 or 4.7 and 5 for this. Target just v5 means not having to worry about making it nice in both, but will mean it'll be a bit more of a delay before release).

I'll admit I've not kept up as much with the 5.0 changes as I should have, but I'm happy either way; I can definitely see the benefit to targeting 5.0 and keeping code maintenance down, or to getting the feature out for folks to use sooner at the cost of more development work.

• I'd like to think about removing Seq and Par and replacing them with some simple rules for sequences, strings and numbers.

Honestly, agreed; I tried a few different designs to try and get rid of those two classes, but they all felt a bit awkward and special-cased. Happy to revise, though, to lower the barrier to using the new feature.

This would match, e.g., qutip.ket("01") and qutip.basis([2, 2], [0, 1]). Seq and Par do however make it really clear that in one case measurements follow each other on the same subspaces and in the other they are performed simultaneously on different subspaces, so I'm not quite sure. Maybe there is some middle ground.

My initial thought was to do something like use tuples instead of Seq and lists instead of Par, but that then ran into a couple issues. Namely, it wasn't as obvious what each kind of container meant, and it was harder to automatically flatten (e.g.: Seq(1, Seq(2, 3), 4) is identical to Seq(1, 2, 3, 4) and Seq(1, Par(2), 3) is identical to Seq(1, 2, 3) but distinct from Seq(1, Par(2, 3))).

One alternative may be to have a single subclass Outcome of tuple instead of both, then overload operators like | to mean parallel. Internally, could use two subclasses like Seq and Par but display visually as Outcome(1, 2 | 3, 4) instead of Outcome(1, Par(2, 3), 4).

• In QuTiP it's more normal to have the helper constructors not on the class as classmethods. I do like secondary constructors as classmethods, but I think in cases like this where there are essentially an infinite number of possible constructors, it makes sense to not "bless" any of them by sticking them on the class.

Sounds good, will move those out to be ordinary functions, then.

• It would be good to have one really nice use case that we could turn into a guide entry in the documentation. The current small examples are great, but it would be good to add one slightly bigger worked example that did something more "exciting".

I've played around with using if_ to compute the noisy channel one gets for state teleportation when operations are perfect but classical outcomes get scrambled; would that be useful, perhaps?

• I would not push everything into the qutip namespace (largely because we would likely not want that in v5, although I should actually check what v5 is doing in qutip.__init__ these days).

Fair enough, easy to drop that as well. Would you want everything dropped, or would some things like the core type QInstrument itself be reasonable to keep in qutip.__all__?

Hoping to hear comments from others!

Thanks for your review, I'll work on addressing, then excited to hear what others think as well! 💕

[hodgestar]

I've played around with using if_ to compute the noisy channel one gets for state teleportation when operations are perfect but classical outcomes get scrambled; would that be useful, perhaps?

That sounds like a great example!

Fair enough, easy to drop that as well. Would you want everything dropped, or would some things like the core type QInstrument itself be reasonable to keep in qutip.all

Let's keep QInstrument in since that feels at the same level as Qobj and QobjEvo and we can think about other things on a case by case basis (and see what v5 is doing).

[cgranade]

@cgranade When we initially talked about this, we also spoke about storing Kraus superoperators with it, but it's not completely clear to me whether this route covers that case. "No" is fine (we don't have to make this work for everything) but if the answer is "Yes" or "Maybe", what would that look like?

At the moment, no. Following our discussion, it sounded like separating the two features may make the most sense such that I wanted to focus first on representing instruments.

I've played around with using if_ to compute the noisy channel one gets for state teleportation when operations are perfect but classical outcomes get scrambled; would that be useful, perhaps?

That sounds like a great example!

Awesome, I'll go add that as a draft PR to the notebooks repo, then, so as to develop both in parallel.

Fair enough, easy to drop that as well. Would you want everything dropped, or would some things like the core type QInstrument itself be reasonable to keep in qutip.all

Let's keep QInstrument in since that feels at the same level as Qobj and QobjEvo and we can think about other things on a case by case basis (and see what v5 is doing).

Sounds good, will do!

[cgranade]

@hodgestar Apologies for taking as long, but went on and rebased my PR to 5.0 and addressed some of your comments. In particular, factory methods have been moved out into ordinary functions, and I added a simple string-based format for outcome labels (leaving Seq and Par for more complex cases as needed).

For the notebook, I ran into the slight issue that the qutip-notebooks repo uses a more copyleft license, but in the meantime I put together a notebook at https://gist.github.com/cgranade/7c2a5a0827dddc4281666ad45763b1ec that includes a few examples of the API in use, including for modeling a simple teleportation channel.

I still need to address PEP8 issues and add tests, but I think it should be a bit further along; thanks for all your help and feedback! 💕

[hodgestar]

@cgranade The new dimensions support has now been merged, so I think we can pick up this PR again. I also did all the running around necessary to re-license the notebooks (qutip/qutip-notebooks#146). Note that the home for new notebooks is https://github.com/qutip/qutip-tutorials.