Merge pull request #136 from MSRudolph/dev

Performance fixes:
- dispatch to Base.iterate() for PauliSum
- replace Pauli rotation gate with gate bitmask in commutes()
- inline iterate(), add!(), sett!()

Author: MSRudolph

Project.toml

@@ -1,7 +1,7 @@
 name = "PauliPropagation"
 uuid = "293282d5-3c99-4fb6-92d0-fd3280a19750"
 authors = ["Manuel S. Rudolph"]
-version = "0.7.0"
+version = "0.7.1"
 
 [deps]
 AcceleratedKernels = "6a4ca0a5-0e36-4168-a932-d9be78d558f1"

README.md

@@ -31,18 +31,6 @@ Pkg.add(url="https://github.com/MSRudolph/PauliPropagation.jl.git", rev="branchn
 where you can use the keyword `rev="branchname"` to install development versions of the package.
 We don't recommend using branches other than `main` or `dev`.
 
-### Clone the repository and install locally 
-Navigate to a local directory where you want to clone this repository into and run the following in a terminal
-```bash
-git clone git@github.com:MSRudolph/PauliPropagation.jl.git
-```
-Inside this cloned repository, you can now freely import `PauliPropagation` or install it into your environment.\
-Alternatively, you can push the relative path to the cloned repository to the Julia package load path called `LOAD_PATH` via
-```julia
-rel_path = "your/relative/path/PauliPropagation"
-push!(LOAD_PATH,rel_path);
-```
-This may require that you have no global installation of `PauliPropagation` in your enviroment.
 
 ### A note on installing Julia 
 It is recommended to install julia using `juliaup` with instructions from [here](https://github.com/JuliaLang/juliaup). Then, Julia's _long-term support_ version (currently a `1.10` version) can be installed via
@@ -139,20 +127,19 @@ Therefore, the trace is equivalent to the sum over the coefficients of Pauli str
 ```
 
 ## Important Notes and Caveats
-- Circuits are specified in the _Schrödinger_ picture, as if operated upon states. Behind the scenes, `propagate()` will (by default) apply the _adjoint_ circuit upon the passed `PauliSum` which is treated as the observable operator.
-- Schrödinger propagation is planned but not yet supported _except_ through manually passing the _adjoint_ of the intended circuit to `propagate()`. This is often easy. For instance, with the circuit order reversed, angles in `PauliRotation` gates are negated, and `CliffordGate` are passed to `transposecliffordmap()`.
+- Circuits are specified in the _Schrödinger_ picture, as if operated upon states. Behind the scenes, `propagate()` will (by default) apply the _adjoint_ circuit upon the passed `PauliSum` which is treated as the observable operator. The default can be changed by passing `heisenberg=false` to `propagate()`, though it will not make simulating dense quantum states efficient. 
+- Schrödinger propagation via `heisenberg=false` is supported since version `0.7`, but not for all gates. So far, we natively support `PauliRotation`, `CliffordGate`, and `<:PauliNoise` gates. `ImaginaryPauliRotation` is _only_ supported with `heisenberg=false`.
 - While Pauli propagation can, in principle, be used for _extended_ stabilizer simulation, we do not currently support sub-exponential strong simulation of stabilizer states.
-- Sampling quantum states is currently not supported.
+- Sampling quantum states is currently not supported, but is coming soon.
 - Many underlying data structures and functions can be used for other purposes involving Pauli operators.
 
 All of the above can be addressed by writing the additional missing code due to the nice extensibility of Julia.
 
 ## Upcoming Features
 This package is still work-in-progress. You will probably find certain features that you would like to have and that are currently missing.\
 Here are some features that we want to implement in the future. Feel free to contribute!
-- **Multi-threading and improved scalability**. Currently, PauliPropagation.jl works uses a single CPU thread and may run your hardware out of memory. Future versions should be even faster and with options to trade-off computational runtime and memory requirements. 
-- **Easier Schrödinger picture propagation**. Currently, the default is Heisenberg and there is no easy way to transpose the gates.
-- **A fast and flexible Surrogate version**. Currently, we provide a version of the Pauli propagation Surrogate that is _good_ and _works_, at least for Pauli gates and Clifford gates. Stay tuned for a whole lot more.
+- **GPU acceleration**. Since version `0.7`, we provide a PauliPropagationCUDA extension in `ext/`. So far, it only works with `PauliRotation` gates and is not yet maximally performant. 
+- **Stochastic evolution**. Propagation methods are mainly memory-limited. We aim to change this and introduce time vs memory trade-offs.
 
 ## How to contribute
 We have a Slack channel `#pauli-propagation` in the [Julia Slack](https://join.slack.com/t/julialang/shared_invite/zt-2zljxdwnl-kSXbwuwFHeERyxSD3iFJdQ).
@@ -187,13 +174,15 @@ If you are publishing research using `PauliPropagation.jl`, please cite this lib
 ```
 
 ## Related publications
-Some of the developers of this package are co-authors in the following papers using Pauli propagation and (at least parts of) this code:
+Some of the developers of this package are co-authors in the following papers using Pauli propagation and (at least parts of) this code.
+If you are using our package, please consider citing some of these works:
 - [Classical simulations of noisy variational quantum circuits](https://arxiv.org/abs/2306.05400)
 - [Classical surrogate simulation of quantum systems with LOWESA](https://arxiv.org/abs/2308.09109)
 - [Quantum Convolutional Neural Networks are (Effectively) Classically Simulable](https://arxiv.org/abs/2408.12739)
 - [Classically estimating observables of noiseless quantum circuits](https://arxiv.org/abs/2409.01706)
 - [Efficient quantum-enhanced classical simulation for patches of quantum landscapes](https://arxiv.org/abs/2411.19896)
 - [Simulating quantum circuits with arbitrary local noise using Pauli Propagation](https://arxiv.org/abs/2501.13101)
 - [Leveraging Symmetry Merging in Pauli Propagation](https://arxiv.org/abs/2512.12094)
+- [Thermal State Simulation with Pauli and Majorana Propagation](https://www.arxiv.org/abs/2602.04878)
 
 And more are coming up.

src/Base/termsum.jl

@@ -95,7 +95,17 @@ end
 
 
 ### Default functions defined for all TermSum types
-function Base.iterate(term_sum::AbstractTermSum)
+@inline function Base.iterate(term_sum::AbstractTermSum, args...)
+    return _iterate(StorageType(term_sum), term_sum, args...)
+end
+
+@inline function _iterate(::DictStorage, term_sum::AbstractTermSum, args...)
+    dict_storage = storage(term_sum)
+    return iterate(dict_storage, args...)
+end
+
+
+@inline function _iterate(::StorageType, term_sum::AbstractTermSum)
     # 1. Create the iterator we are delegating to
     iter = zip(terms(term_sum), coefficients(term_sum))
 
@@ -107,7 +117,7 @@ function Base.iterate(term_sum::AbstractTermSum)
     return next === nothing ? nothing : (next[1], (iter, next[2]))
 end
 
-function Base.iterate(term_sum::AbstractTermSum, state)
+@inline function _iterate(::StorageType, term_sum::AbstractTermSum, state)
     # 1. Unpack the state tuple
     (iter, inner_state) = state
 
@@ -138,16 +148,16 @@ end
     add!(term_sum::AbstractTermSum, term, coeff)
 
 Add `coeff` to the coefficient of `term` in `term_sum`.
-Calls `add!(StorageType(term_sum), term_sum, term, coeff)` internally.
-For custom behavior, overload `storage()` and/or `add!` for the specific TermSum type.
+Calls `_add!(StorageType(term_sum), term_sum, term, coeff)` internally.
+For custom behavior, overload `storage()` and/or `_add!` for the specific TermSum type.
 """
-function add!(term_sum::AbstractTermSum, term, coeff)
-    add!(StorageType(term_sum), term_sum, term, coeff)
+@inline function add!(term_sum::AbstractTermSum, term, coeff)
+    _add!(StorageType(term_sum), term_sum, term, coeff)
     return term_sum
 end
 
 
-function add!(::DictStorage, term_sum::AbstractTermSum, term, coeff)
+@inline function _add!(::DictStorage, term_sum::AbstractTermSum, term, coeff)
     dict_storage = storage(term_sum)
     if haskey(dict_storage, term)
         dict_storage[term] += coeff
@@ -157,7 +167,7 @@ function add!(::DictStorage, term_sum::AbstractTermSum, term, coeff)
     return term_sum
 end
 
-function add!(::ArrayStorage, term_sum::AbstractTermSum, term, coeff)
+@inline function _add!(::ArrayStorage, term_sum::AbstractTermSum, term, coeff)
     terms_vec, coeffs_vec = storage(term_sum)
     ind = findfirst(t -> t == term, terms_vec)
     if !isnothing(ind)
@@ -169,20 +179,20 @@ function add!(::ArrayStorage, term_sum::AbstractTermSum, term, coeff)
     return term_sum
 end
 
-function add!(::StorageType, term_sum::AbstractTermSum, term, coeff)
-    thrownotimplemented(typeof(term_sum), :add!)
+@inline function _add!(::StorageType, term_sum::AbstractTermSum, term, coeff)
+    thrownotimplemented(typeof(term_sum), :_add!)
 end
 
 
 function add!(term_sum1::AbstractTermSum, term_sum2::AbstractTermSum)
-    for (term, coeff) in zip(terms(term_sum2), coefficients(term_sum2))
+    for (term, coeff) in term_sum2
         add!(term_sum1, term, coeff)
     end
     return term_sum1
 end
 
 function subtract!(term_sum1::AbstractTermSum, term_sum2::AbstractTermSum)
-    for (term, coeff) in zip(terms(term_sum2), coefficients(term_sum2))
+    for (term, coeff) in term_sum2
         add!(term_sum1, term, -coeff)
     end
     return term_sum1
@@ -193,23 +203,22 @@ end
     set!(term_sum::AbstractTermSum, term, coeff)
     
 Set the coefficient of `term` in `term_sum` to `coeff`.
-Calls `set!(StorageType(term_sum), term_sum, term, coeff)` internally.
-For custom behavior, overload `storage()` and/or `set!` for the specific TermSum type.
+Calls `_set!(StorageType(term_sum), term_sum, term, coeff)` internally.
+For custom behavior, overload `storage()` and/or `_set!` for the specific TermSum type.
 """
-set!(term_sum::AbstractTermSum, term, coeff) = set!(StorageType(term_sum), term_sum, term, coeff)
+@inline function set!(term_sum::AbstractTermSum, term, coeff)
+    _set!(StorageType(term_sum), term_sum, term, coeff)
+    return term_sum
+end
 
-"""
-    set!(term_storage::Dict, term, coeff)
 
-Default implementation of `set!()` for `Dict` storage.
-"""
-function set!(::DictStorage, term_sum::AbstractTermSum, term, coeff)
+@inline function _set!(::DictStorage, term_sum::AbstractTermSum, term, coeff)
     dict_storage = storage(term_sum)
     dict_storage[term] = coeff
     return term_sum
 end
 
-function set!(::ArrayStorage, term_sum::AbstractTermSum, term, coeff)
+@inline function _set!(::ArrayStorage, term_sum::AbstractTermSum, term, coeff)
     terms_vec, coeffs_vec = storage(term_sum)
     ind = findfirst(t -> t == term, terms_vec)
     if !isnothing(ind)
@@ -221,13 +230,13 @@ function set!(::ArrayStorage, term_sum::AbstractTermSum, term, coeff)
     return term_sum
 end
 
-function set!(::StorageType, term_sum::AbstractTermSum, term, coeff)
+@inline function _set!(ST::StorageType, term_sum::AbstractTermSum, term, coeff)
     old_coeff = getcoeff(term_sum, term)
     if old_coeff == zero(coefftype(term_sum))
-        add!(term_sum, term, coeff)
+        _add!(ST, term_sum, term, coeff)
     else
         delta = coeff - old_coeff
-        add!(term_sum, term, delta)
+        _add!(ST, term_sum, term, delta)
     end
     return term_sum
 end

src/Propagation/specializations.jl

@@ -30,7 +30,7 @@ function PropagationBase.applytoall!(gate::PauliRotation, prop_cache::PauliPropa
     # loop over all Pauli strings and their coefficients in the Pauli sum
     for (pstr, coeff) in psum
 
-        if commutes(gate, pstr)
+        if commutes(gate_mask, pstr)
             # if the gate commutes with the pauli string, do nothing
             continue
         end

src/Visualization/treetracker.jl

@@ -175,12 +175,12 @@ function PropagationBase.applytoall!(gate::PauliRotation, prop_cache::PauliPropa
     aux_psum = auxsum(prop_cache)
 
     # Convert PauliRotation to MaskedPauliRotation for efficiency
-    gate_mask = symboltoint(nqubits(psum), gate.symbols, gate.qinds)
+    gate_mask = symboltoint(paulitype(psum), gate.symbols, gate.qinds)
 
     # Loop over all Pauli strings and their coefficients in the Pauli sum
     for (pstr, coeff) in psum
         pstr_str = format_pauli_string(pstr, psum.nqubits)
-        if commutes(gate, pstr)
+        if commutes(gate_mask, pstr)
             # If the gate commutes, create a new child node and edge
             pstr_str = format_pauli_string(pstr, psum.nqubits)
             gate_symbol = isempty(gate.symbols) ? "?" : string(gate.symbols[1])