.. _bosonic_synthesis_explanation:

Synthesize an Evolution gate in bosonic space
=============================================

.. meta::
   :description: How to synthesize an Evolution gate into smaller bosonic gates with product
                 formulas (Trotter and qDRIFT) before any truncation or encoding is chosen, on one
                 gate or across a whole circuit, and why this reduces the error between
                 non-commuting groups.

The time evolution :math:`e^{-itH}` of a Hamiltonian that is a sum of non-commuting pieces
cannot be implemented exactly gate by gate. A **product formula** approximates it by a sequence
of short evolutions under the individual pieces, and the approximation error is called the
Trotter error.

You can apply a product formula in two places, and they are not interchangeable:

1. **In bosonic space**, before any truncation or encoding has been chosen. The formula sees the
   *groups* of the operator (see :ref:`grouping_explanation`) and can interleave them.
2. **In qubit space**, inside Qiskit's
   :class:`~qiskit.circuit.library.PauliEvolutionGate`, after the operator has been mapped. The
   formula sees one Pauli operator per group and can only subdivide the terms *within* it.

Only the first can reduce the error between groups that do not commute, which is usually the
dominant one. This guide shows how to use it.

.. important::
   Both steps are optional. An :class:`~qiskit_bosons.circuit.library.Evolution` gate that is
   never decomposed in bosonic space reaches the boson-to-qubit stage whole and is mapped and
   synthesized there. That incurs *no* error at the bosonic level, only whatever the qubit-level
   formula introduces.

Set up
------

A two-site Bose-Hubbard chain, the smallest example in which the groups do not commute (see
:ref:`bose_hubbard_getting_started` for the full model). It is grouped into one Hermitian group for
the bond and one diagonal group for the on-site interaction:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> import numpy as np
   >>> from qiskit.quantum_info import SparsePauliOp, Statevector
   >>> from qiskit.transpiler import PassManager
   >>> from scipy.sparse import csc_matrix
   >>> from scipy.sparse.linalg import expm_multiply
   >>> from qiskit_bosons.circuit import BosonicCircuit
   >>> from qiskit_bosons.circuit.library import Evolution, InitializeModes
   >>> from qiskit_bosons.mappers import LevelSpec
   >>> from qiskit_bosons.mappers.library import occupation_encoding
   >>> from qiskit_bosons.operators import BosonOperator, ann, cre
   >>> from qiskit_bosons.transpiler.presets import generate_preset_occupation_pass_manager
   >>>
   >>> num_sites, max_occupation, total_time = 2, 2, 0.8
   >>> bonds = [(0, 1)]
   >>> terms = {}
   >>> for left, right in bonds:
   ...     terms[(cre(left), ann(right))] = -1.0
   ...     terms[(cre(right), ann(left))] = -1.0
   >>> for site in range(num_sites):
   ...     terms[(cre(site), cre(site), ann(site), ann(site))] = 1.0
   >>> hamiltonian = BosonOperator.from_dict(terms)
   >>> hamiltonian.groups = [
   ...     len(bonds)
   ...     if len({mode for _, mode in term}) == 1
   ...     else bonds.index(tuple(sorted({mode for _, mode in term})))
   ...     for term, _ in hamiltonian.iter_terms()
   ... ]
   >>> hamiltonian.num_groups()
   2

To measure an error you need a reference. The exact evolution of the Hamiltonian *as mapped at the
same truncation* is the right one, because the truncation is a separate approximation:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> spec = LevelSpec(max_occupation, "unary_projected")
   >>> mapped = occupation_encoding(hamiltonian, spec).simplify()
   >>> matrix = csc_matrix(SparsePauliOp.from_sparse_observable(mapped).to_matrix(sparse=True))
   >>> pm = generate_preset_occupation_pass_manager(spec)
   >>> pm.qubit = PassManager()  # keep each group as one exact PauliEvolutionGate
   >>>
   >>> def build(**evolution_kwargs):
   ...     circuit = BosonicCircuit(num_sites)
   ...     circuit.append(InitializeModes([1, 1]), circuit.modes)
   ...     circuit.append(
   ...         Evolution(num_sites, hamiltonian, time=total_time, **evolution_kwargs), circuit.modes
   ...     )
   ...     return circuit
   >>>
   >>> # Under the unary encodings, occupation ``n`` of mode ``m`` sets qubit ``m * d + n``, with
   >>> # ``d = max_occupation + 1`` levels per mode. That places |1, 1> on the basis state below.
   >>> levels = max_occupation + 1
   >>> index = sum(1 << (mode * levels + n) for mode, n in enumerate([1, 1]))
   >>> initial = Statevector.from_int(index, 2**(num_sites * levels)).data
   >>> exact = expm_multiply(-1j * total_time * matrix, initial)

   >>> def error(circuit):
   ...     """Distance to the exact evolution, of the circuit's product of group exponentials."""
   ...     lowered = pm.run(circuit)
   ...     state = Statevector.from_int(0, 2**lowered.num_qubits).evolve(lowered).data
   ...     return float(np.max(np.abs(state - exact)))

.. note::
   ``error`` deliberately simulates the lowered circuit *without* decomposing its
   :class:`~qiskit.circuit.library.PauliEvolutionGate` instances.
   :class:`~qiskit.quantum_info.Statevector` applies each as the exact exponential of its group, so the only thing left is the product of group exponentials that the bosonic synthesis chose.

The default
-----------

With no ``synthesis`` argument an :class:`~qiskit_bosons.circuit.library.Evolution` gate
decomposes with :class:`~qiskit_bosons.circuit.library.synthesis.BosonicLieTrotter`, the
first-order formula: one factor per group, each for the whole time.

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> dict(build().decompose().count_ops())
   {'Evolution': 2, 'InitializeModes': 1}

Without decomposing, the whole gate is handed to the qubit stage and there is no bosonic Trotter
error at all:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> error(build()) < 1e-6
   True

Decomposed, a genuine error appears, because the two groups do not commute:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> error(build().decompose()) > 0.1
   True

Reduce the error
----------------

Pass a different synthesis method to the gate. A higher **order** composes the factors
symmetrically so that the leading error terms cancel; more **repetitions** divide the time into
shorter steps.

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_bosons.circuit.library.synthesis import BosonicLieTrotter, BosonicSuzukiTrotter
   >>>
   >>> for label, synthesis in [
   ...     ("Lie-Trotter", BosonicLieTrotter()),
   ...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
   ...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
   ...     ("Suzuki, order 4*", BosonicSuzukiTrotter(order=4, merge_adjacent=False)),
   ... ]:
   ...     circuit = build(synthesis=synthesis).decompose()
   ...     print(f"{label:16s} {circuit.count_ops()['Evolution']:3d} factors")
   Lie-Trotter        2 factors
   Suzuki, order 2    3 factors
   Suzuki, order 4   11 factors
   Suzuki, order 4*  15 factors

Each increase of the order by two composes five sweeps of the one below (the :math:`5^{(k-2)/2}`
growth of Suzuki's recursion, applied to the :math:`2n - 1` factors of order 2), which is the count
in the last line, marked with a star. The sweeps meet on the same factor, so the neighboring gates
evolve under one operator and are fused into a single gate with the summed time. That saves one gate
per seam and is exact. Pass ``merge_adjacent=False`` to keep the textbook sequence. Either way, a
higher order buys a smaller error:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> lie = error(build(synthesis=BosonicLieTrotter()).decompose())
   >>> second = error(build(synthesis=BosonicSuzukiTrotter(order=2)).decompose())
   >>> fourth = error(build(synthesis=BosonicSuzukiTrotter(order=4)).decompose())
   >>> lie > second > fourth
   True

Repetitions trade depth for accuracy at a different rate. The first-order error falls linearly in
the step size, the second-order one quadratically:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> lie_errors = [error(build(synthesis=BosonicLieTrotter(reps=r)).decompose()) for r in (1, 2, 4, 8)]
   >>> bool(all(a > b for a, b in zip(lie_errors, lie_errors[1:])))
   True
   >>> bool(lie_errors[0] / lie_errors[-1] > 6)  # about 8, i.e. linear in 1 / reps
   True
   >>> suzuki_errors = [
   ...     error(build(synthesis=BosonicSuzukiTrotter(order=2, reps=r)).decompose()) for r in (1, 2, 4)
   ... ]
   >>> bool(suzuki_errors[0] / suzuki_errors[-1] > 10)  # about 16, i.e. quadratic in 1 / reps
   True

Which is the better trade depends on the operator. Note that an operator with a single group is
exact at any order, because a lone factor commutes with itself, and the formula then collapses to one
gate rather than re-splitting it.

Measure the cost
----------------

A smaller error costs more factors, and each factor costs entangling gates once the qubit stage has
synthesized it. How many depends on the encoding, so count them through the preset pipeline, with
the basis fixed and no qubit-level optimization (``optimization_level=0``), to see the price of each
method:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> basis = ["rz", "sx", "x", "cx"]
   >>>
   >>> def cx_count(circuit, encoding):
   ...     spec = LevelSpec(max_occupation, encoding)
   ...     lowered = generate_preset_occupation_pass_manager(
   ...         spec, optimization_level=0, basis_gates=basis
   ...     ).run(circuit)
   ...     return lowered.count_ops().get("cx", 0)
   >>>
   >>> methods = [
   ...     ("Lie-Trotter", BosonicLieTrotter()),
   ...     ("Suzuki, order 2", BosonicSuzukiTrotter(order=2)),
   ...     ("Suzuki, order 4", BosonicSuzukiTrotter(order=4)),
   ... ]
   >>> encodings = ("gray", "unary_unprojected")
   >>> counts = {
   ...     encoding: [
   ...         cx_count(build(synthesis=method).decompose(), encoding) for _, method in methods
   ...     ]
   ...     for encoding in encodings
   ... }
   >>> # the order buys accuracy at a price that grows with it
   >>> bool(counts["gray"][0] < counts["gray"][1] < counts["gray"][2])
   True

The exact counts depend on how Qiskit synthesizes the evolution gates, which changes between
releases, so the plot shows them rather than a table:

.. plot::
   :context: close-figs
   :alt: Grouped bar chart of the CX count for Lie-Trotter and Suzuki of order 2 and 4, for the Gray and the unary encoding. The counts grow with the order, and the unary bars are higher.

   >>> import matplotlib.pyplot as plt
   >>>
   >>> figure, axes = plt.subplots(figsize=(6, 3.5), layout="constrained")
   >>> width = 0.38
   >>> positions = list(range(len(methods)))
   >>> for offset, encoding in zip((-width / 2, width / 2), encodings):
   ...     _ = axes.bar([p + offset for p in positions], counts[encoding], width, label=encoding)
   >>> _ = axes.set_xticks(positions, [name for name, _ in methods])
   >>> _ = axes.set(ylabel="CX count")
   >>> _ = axes.legend()
   >>> figure
   <Figure size ... with 1 Axes>

Each group of bars is a method: Lie-Trotter, Suzuki of order 2, and Suzuki of order 4. The order buys
accuracy at a price that grows by the factor of the recursion, and the unary encoding pays more than
the Gray code for every method because each of its factors maps to more entangling gates. Choose the
method for the accuracy that you need and the encoding for the hardware, and keep the two choices
apart: the bosonic circuit above never changes.

Apply one method to a whole circuit
-----------------------------------

To set ``synthesis`` on each gate, you must thread the choice through everything that builds one.
The :class:`~qiskit_bosons.transpiler.passes.B2BSynthesis` pass makes the choice once for a
pipeline instead. Its ``config`` names the method for each type of gate, and it both selects the
method and expands the gates:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> from qiskit_bosons.transpiler import BosonicPassManager
   >>> from qiskit_bosons.transpiler.passes import B2BSynthesis
   >>>
   >>> trotterize = BosonicPassManager([B2BSynthesis({"Evolution": BosonicSuzukiTrotter(order=2)})])
   >>> dict(trotterize.run(build()).count_ops())
   {'Evolution': 3, 'InitializeModes': 1}

A method can also be named, with the arguments for its ``__init__``, instead of being built:
``B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})})``. The available names are those
of the :class:`~qiskit_bosons.circuit.library.synthesis.B2BSynthesisPluginManager`, and a method that
is registered nowhere can be handed over as an instance, as above.

With ``apply=False`` the pass only selects the method and leaves the expansion to the point where
it is wanted, for example to inspect the circuit first. An atomic gate is never touched, not even
tagged:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> selected = BosonicPassManager(
   ...     [B2BSynthesis({"Evolution": ("SuzukiTrotter", (), {"order": 2})}, apply=False)]
   ... ).run(build())
   >>> selected.count_ops()["Evolution"]  # selected, not yet expanded
   1
   >>> selected.decompose().count_ops()["Evolution"]  # expanded on demand
   3

An unexpanded gate reaches the boson-to-qubit stage whole, where it is mapped without its
``synthesis`` being consulted, so every method then produces identical output.

An optional ``filter`` restricts it to the gates that benefit. Evolutions whose groups all
commute, such as a Hamiltonian made only of number operators, are synthesized exactly at any
order, so a higher order would only add depth there.

Restrict the pass to the gates that benefit
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

A purely diagonal operator makes the point concrete, because number operators commute with each
other. Its groups are exact at first order, and the second order spends entangling gates for
nothing:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> # one group per term, each Hermitian; the order is fixed because order two evolves the first
   >>> # group twice, so which group comes first changes the cost
   >>> diagonal = BosonOperator.from_terms_with_groups(
   ...     [([cre(0), ann(0), cre(1), ann(1)], 1.0, 0), ([cre(0), ann(0)], 0.5, 1)]
   ... )
   >>>
   >>> def diagonal_cost(synthesis):
   ...     circuit = BosonicCircuit(num_sites)
   ...     circuit.append(Evolution(num_sites, diagonal, time=total_time, synthesis=synthesis), circuit.modes)
   ...     return cx_count(circuit.decompose(), "gray")
   >>>
   >>> bool(diagonal_cost(BosonicSuzukiTrotter(order=2)) > diagonal_cost(BosonicLieTrotter()))
   True

Both circuits are exact, so the second spends more entangling gates for no gain. (The
symmetric order-two formula applies the first group twice at half the time each, and thus the
cost depends on which group comes first.) The
``filter`` takes a predicate over the :class:`~qiskit.dagcircuit.DAGOpNode` and leaves the rejected
gates untouched, so one pass can raise the order on the gates that benefit and keep the rest at
first order. The predicate below accepts a gate only when its operator moves bosons between modes,
which is when the groups fail to commute:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> def is_diagonal(operator):
   ...     """Whether every term only counts bosons, rather than moving them between modes."""
   ...     return all(
   ...         sorted(mode for created, mode in actions if created)
   ...         == sorted(mode for created, mode in actions if not created)
   ...         for actions, _ in operator.iter_terms()
   ...     )
   >>>
   >>> def mixed_circuit():
   ...     circuit = BosonicCircuit(num_sites)
   ...     circuit.append(Evolution(num_sites, hamiltonian, time=total_time), circuit.modes)
   ...     circuit.append(Evolution(num_sites, diagonal, time=total_time), circuit.modes)
   ...     return circuit
   >>>
   >>> order_2 = BosonicSuzukiTrotter(order=2)
   >>> everywhere = cx_count(
   ...     BosonicPassManager([B2BSynthesis({"Evolution": order_2})]).run(mixed_circuit()), "gray"
   ... )
   >>> restricted = cx_count(
   ...     BosonicPassManager(
   ...         [
   ...             B2BSynthesis(
   ...                 {"Evolution": order_2},
   ...                 filter=lambda node: not is_diagonal(node.op.operator),
   ...             )
   ...         ]
   ...     ).run(mixed_circuit()),
   ...     "gray",
   ... )
   >>> bool(everywhere > restricted)  # the filter saves entangling gates
   True

If you run the pass again, nothing changes: the factors it emits are *atomic*.

Atomic factors
--------------

An :class:`~qiskit_bosons.circuit.library.Evolution` marked ``atomic`` has no definition, so
:meth:`~qiskit_bosons.circuit.BosonicCircuit.decompose` leaves it in place however often it is
applied. A synthesis marks the factors it emits this way, and so does
:class:`~qiskit_bosons.circuit.library.synthesis.BosonicQDrift` for the terms it samples.

This is more than tidiness. Splitting a Hermitian group further, term by term, yields factors
that are individually **not** Hermitian, and the exponential of a non-Hermitian operator is not
unitary:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> circuit = build()
   >>> once = circuit.decompose()
   >>> dict(circuit.decompose(reps=5).count_ops()) == dict(once.count_ops())
   True

.. caution::
   Each factor must be Hermitian on its own. A Hermitian operator split into groups can produce
   non-Hermitian ones, such as a hopping term separated from its adjoint. Check with
   :func:`~qiskit_bosons.operators.terms.groups_are_hermitian`. An operator *without* groups is
   split term by term, and the qubit stage then rejects the factors with ``Operator contains
   complex coefficients``. The symmetrization of an even order partially cancels the error of a
   non-Hermitian factor, so a wrongly grouped operator can look better at order 2 than at order 1
   while still being wrong.

Repeat a short step
-------------------

For a circuit-level alternative, build one short Trotter step and repeat it.
:meth:`~qiskit_bosons.circuit.BosonicCircuit.repeat` flattens the repetitions into one circuit and
leaves the boundaries adjacent, so passes can act across them:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> step = BosonicCircuit(num_sites)
   >>> step.append(Evolution(num_sites, hamiltonian, time=total_time / 4), step.modes)
   >>> dict(step.repeat(4).count_ops())
   {'Evolution': 4}

Pass ``insert_barriers=True`` to keep the repetitions apart. A barrier is mapped onto every qubit
of the registers it touches, because an encoding needs several qubits per mode and a mode has no
per-qubit image of its own:

.. plot::
   :context:
   :nofigs:
   :include-source:

   >>> [inst.operation.name for inst in step.repeat(3, insert_barriers=True)._inner.data]
   ['Evolution', 'barrier', 'Evolution', 'barrier', 'Evolution']

.. note::
   ``repeat`` is unrelated to the ``reps`` of a synthesis method. A synthesis subdivides a
   *single* evolution into more, shorter factors. Repeating a circuit that holds one evolution
   merely multiplies its total time.

Combine with a qubit-level formula
----------------------------------

The two levels compose, and where both approximate, the accuracy of the result is governed by the
weaker of the two. Raising the bosonic order while the qubit-side formula stays first-order buys
little. See :class:`~qiskit_bosons.transpiler.passes.MapperFnEvolutionSynthesis` for the
qubit-side formula. It keeps the Pauli term order that the mapper chose, and so
:func:`~qiskit_bosons.transpiler.passes.group_wise` can emit the terms of each group adjacent to one
another, and :func:`~qiskit_bosons.transpiler.passes.simplify` can pin a canonical order instead.

Next steps
----------

- :ref:`transpilation_explanation` explains the stages that follow the bosonic synthesis, and
  where :class:`~qiskit_bosons.transpiler.passes.B2BSynthesis` sits among them.
- :ref:`grouping_explanation` explains how to assign and check the groups that a product formula
  interleaves.
- :ref:`bose_hubbard_getting_started` applies the repetition scheme to a complete example.