Migrating the Python operators to be C++ bindings

docs/sphinx/examples/python/dynamics/cavity_qed.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import operators, Schedule, ScipyZvodeIntegrator
+from cudaq import boson, Schedule, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os
@@ -15,17 +15,17 @@
 
 # Alias for commonly used operators
 # Cavity operators
-a = operators.annihilate(1)
-a_dag = operators.create(1)
+a = boson.annihilate(1)
+a_dag = boson.create(1)
 
 # Atom operators
-sm = operators.annihilate(0)
-sm_dag = operators.create(0)
+sm = boson.annihilate(0)
+sm_dag = boson.create(0)
 
 # Defining the Hamiltonian for the system: self-energy terms and cavity-atom interaction term.
 # This is the so-called Jaynes-Cummings model:
 # https://en.wikipedia.org/wiki/Jaynes%E2%80%93Cummings_model
-hamiltonian = 2 * np.pi * operators.number(1) + 2 * np.pi * operators.number(
+hamiltonian = 2 * np.pi * boson.number(1) + 2 * np.pi * boson.number(
     0) + 2 * np.pi * 0.25 * (sm * a_dag + sm_dag * a)
 
 # Initial state of the system
@@ -41,23 +41,23 @@
 schedule = Schedule(steps, ["time"])
 
 # First, evolve the system without any collapse operators (ideal).
-evolution_result = cudaq.evolve(
-    hamiltonian,
-    dimensions,
-    schedule,
-    rho0,
-    observables=[operators.number(1), operators.number(0)],
-    collapse_operators=[],
-    store_intermediate_results=True,
-    integrator=ScipyZvodeIntegrator())
+evolution_result = cudaq.evolve(hamiltonian,
+                                dimensions,
+                                schedule,
+                                rho0,
+                                observables=[boson.number(1),
+                                             boson.number(0)],
+                                collapse_operators=[],
+                                store_intermediate_results=True,
+                                integrator=ScipyZvodeIntegrator())
 
 # Then, evolve the system with a collapse operator modeling cavity decay (leaking photons)
 evolution_result_decay = cudaq.evolve(
     hamiltonian,
     dimensions,
     schedule,
     rho0,
-    observables=[operators.number(1), operators.number(0)],
+    observables=[boson.number(1), boson.number(0)],
     collapse_operators=[np.sqrt(0.1) * a],
     store_intermediate_results=True,
     integrator=ScipyZvodeIntegrator())

docs/sphinx/examples/python/dynamics/gate_calibration.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import operators, Schedule, ScalarOperator, ScipyZvodeIntegrator
+from cudaq import boson, Schedule, ScalarOperator, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os
@@ -55,16 +55,16 @@ def cost_function(amps):
     amplitude = 100 * amps[0]
     drag_amp = 100 * amps[1]
     # Qubit Hamiltonian
-    hamiltonian = detuning * operators.number(0) + (
-        anharmonicity / 2) * operators.create(0) * operators.create(
-            0) * operators.annihilate(0) * operators.annihilate(0)
+    hamiltonian = detuning * boson.number(0) + (
+        anharmonicity / 2) * boson.create(0) * boson.create(
+            0) * boson.annihilate(0) * boson.annihilate(0)
     # Drive term
-    hamiltonian += amplitude * ScalarOperator(gaussian) * (
-        operators.create(0) + operators.annihilate(0))
+    hamiltonian += amplitude * ScalarOperator(gaussian) * (boson.create(0) +
+                                                           boson.annihilate(0))
 
     # Drag term (leakage reduction)
     hamiltonian += 1j * drag_amp * ScalarOperator(dgaussian) * (
-        operators.annihilate(0) - operators.create(0))
+        boson.annihilate(0) - boson.create(0))
 
     # We optimize for a X(pi/2) rotation
     evolution_result = cudaq.evolve(hamiltonian,

docs/sphinx/examples/python/dynamics/heisenberg_model.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import operators, spin, Schedule, ScipyZvodeIntegrator
+from cudaq import spin, Schedule, ScipyZvodeIntegrator
 
 import numpy as np
 import cupy as cp
@@ -25,7 +25,7 @@
     spin_state += str(int(i % 2))
 
 # Observable is the staggered magnetization operator
-staggered_magnetization_op = operators.zero()
+staggered_magnetization_op = spin.empty()
 for i in range(N):
     if i % 2 == 0:
         staggered_magnetization_op += spin.z(i)
@@ -42,7 +42,7 @@
     Jz = g
 
     # Construct the Hamiltonian
-    H = operators.zero()
+    H = spin.empty()
 
     for i in range(N - 1):
         H += Jx * spin.x(i) * spin.x(i + 1)

docs/sphinx/examples/python/dynamics/landau_zener.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import spin, operators, ScalarOperator, Schedule, ScipyZvodeIntegrator
+from cudaq import spin, boson, ScalarOperator, Schedule, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os
@@ -20,8 +20,8 @@
 # Define some shorthand operators
 sx = spin.x(0)
 sz = spin.z(0)
-sm = operators.annihilate(0)
-sm_dag = operators.create(0)
+sm = boson.annihilate(0)
+sm_dag = boson.create(0)
 
 # Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
 dimensions = {0: 2}
@@ -51,7 +51,7 @@
                                 dimensions,
                                 schedule,
                                 psi0,
-                                observables=[operators.number(0)],
+                                observables=[boson.number(0)],
                                 collapse_operators=[],
                                 store_intermediate_results=True,
                                 integrator=ScipyZvodeIntegrator())

docs/sphinx/examples/python/dynamics/mgmn/initial_state.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import operators, spin, Schedule, RungeKuttaIntegrator
+from cudaq import spin, Schedule, RungeKuttaIntegrator
 
 import numpy as np
 import matplotlib.pyplot as plt
@@ -19,14 +19,14 @@
     dimensions[i] = 2
 
 # Observable is the average magnetization operator
-avg_magnetization_op = operators.zero()
+avg_magnetization_op = spin.empty()
 for i in range(N):
     avg_magnetization_op += (spin.z(i) / N)
 
 # Arbitrary coupling constant
 g = 1.0
 # Construct the Hamiltonian
-H = operators.zero()
+H = spin.empty()
 for i in range(N):
     H += 2 * np.pi * spin.x(i)
     H += 2 * np.pi * spin.y(i)
@@ -38,9 +38,9 @@
 schedule = Schedule(steps, ["time"])
 
 # Initial state (expressed as an enum)
-psi0 = cudaq.operator.InitialState.ZERO
+psi0 = cudaq.dynamics.InitialState.ZERO
 # This can also be used to initialize a uniformly-distributed wave-function instead.
-# `psi0 = cudaq.operator.InitialState.UNIFORM`
+# `psi0 = cudaq.dynamics.InitialState.UNIFORM`
 
 # Run the simulation
 evolution_result = cudaq.evolve(H,

docs/sphinx/examples/python/dynamics/mgmn/multi_gpu.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import operators, spin, Schedule, RungeKuttaIntegrator
+from cudaq import spin, Schedule, RungeKuttaIntegrator
 
 import numpy as np
 import cupy as cp
@@ -29,7 +29,7 @@
     spin_state += str(int(i % 2))
 
 # Observable is the staggered magnetization operator
-staggered_magnetization_op = operators.zero()
+staggered_magnetization_op = spin.empty()
 for i in range(N):
     if i % 2 == 0:
         staggered_magnetization_op += spin.z(i)
@@ -46,7 +46,7 @@
     Jz = g
 
     # Construct the Hamiltonian
-    H = operators.zero()
+    H = spin.empty()
 
     for i in range(N - 1):
         H += Jx * spin.x(i) * spin.x(i + 1)

docs/sphinx/examples/python/dynamics/pulse.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import spin, operators, ScalarOperator, Schedule, ScipyZvodeIntegrator
+from cudaq import spin, boson, ScalarOperator, Schedule, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os
@@ -59,7 +59,7 @@ def signal(t):
                                 dimensions,
                                 schedule,
                                 psi0,
-                                observables=[operators.number(0)],
+                                observables=[boson.number(0)],
                                 collapse_operators=[],
                                 store_intermediate_results=True,
                                 integrator=ScipyZvodeIntegrator())

docs/sphinx/examples/python/dynamics/silicon_spin_qubit.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import spin, operators, Schedule, ScalarOperator, ScipyZvodeIntegrator
+from cudaq import spin, boson, Schedule, ScalarOperator, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os
@@ -36,7 +36,7 @@
                                     dimensions,
                                     schedule,
                                     rho0,
-                                    observables=[operators.number(0)],
+                                    observables=[boson.number(0)],
                                     collapse_operators=[],
                                     store_intermediate_results=True,
                                     integrator=ScipyZvodeIntegrator())

docs/sphinx/examples/python/dynamics/tensor_callback.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import ElementaryOperator, spin, operators, Schedule, ScipyZvodeIntegrator
+from cudaq import MatrixOperatorElement, operators, boson, Schedule, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os
@@ -26,10 +26,10 @@ def callback_tensor(t):
 lz_formula_p1 = 1.0 - lz_formula_p0
 
 # Let's define the control term as a callback tensor that acts on 2-level systems
-ElementaryOperator.define("lz_op", [2], callback_tensor)
+operators.define("lz_op", [2], callback_tensor)
 
 # Hamiltonian
-hamiltonian = ElementaryOperator("lz_op", [0])
+hamiltonian = operators.instantiate("lz_op", [0])
 
 # Dimensions of sub-system. We only have a single degree of freedom of dimension 2 (two-level system).
 dimensions = {0: 2}
@@ -46,7 +46,7 @@ def callback_tensor(t):
                                 dimensions,
                                 schedule,
                                 psi0,
-                                observables=[operators.number(0)],
+                                observables=[boson.number(0)],
                                 collapse_operators=[],
                                 store_intermediate_results=True,
                                 integrator=ScipyZvodeIntegrator())

docs/sphinx/examples/python/dynamics/transmon_resonator.py

@@ -1,5 +1,5 @@
 import cudaq
-from cudaq import operators, spin, Schedule, ScipyZvodeIntegrator
+from cudaq import operators, spin, operators, Schedule, ScipyZvodeIntegrator
 import numpy as np
 import cupy as cp
 import os

python/CMakeLists.txt

@@ -78,4 +78,4 @@ if(CUDAQ_BUILD_TESTS)
 endif()
 
 add_subdirectory(runtime/interop)
-add_subdirectory(cudaq/operator)
+add_subdirectory(cudaq/dynamics)

python/cudaq/__init__.py

@@ -84,18 +84,18 @@
 qreg = cudaq_runtime.qvector
 
 # Operator API
-from .operator.definitions import spin, SpinOperator
-# Re-export non-spin-op operators under the `operators` namespace
-# e.g. `cudaq.operators.annihilate`
-from .operator import operators as operators
-# Operator types (in addition to `spin` types)
-# e.g., allows users to use `cudaq.ScalarOperator(lambda...)`
-from .operator import Operator, ElementaryOperator, ScalarOperator
+from .operators import boson
+from .operators import fermion
+from .operators import spin
+from .operators import custom as operators
+from .operators.definitions import *
+from .operators.manipulation import OperatorArithmetics
+import cudaq.operators.expressions  # needs to be imported, since otherwise e.g. evaluate is not defined
 
 # Time evolution API
-from .operator.schedule import Schedule
-from .operator.evolution import evolve, evolve_async
-from .operator.integrators import *
+from .dynamics.schedule import Schedule
+from .dynamics.evolution import evolve, evolve_async
+from .dynamics.integrators import *
 
 # Optimizers + Gradients
 optimizers = cudaq_runtime.optimizers

python/cudaq/operator/__init__.py → python/cudaq/dynamics/__init__.py

@@ -6,8 +6,5 @@
 # the terms of the Apache License 2.0 which accompanies this distribution.     #
 # ============================================================================ #
 
-from .definitions import operators, spin
-from .evolution import evolve, evolve_async
-from .expressions import Operator, OperatorSum, ProductOperator, ElementaryOperator, ScalarOperator, RydbergHamiltonian
-from .helpers import NumericType, InitialState
-from .schedule import Schedule
+from .helpers import InitialState
+from .schedule import Schedule

python/cudaq/operator/cudm_helpers.py → python/cudaq/dynamics/cudm_helpers.py

@@ -9,10 +9,8 @@
 import logging
 import numpy
 from numbers import Number
-from typing import Any, Mapping, List, Sequence, Union
-
-from .expressions import ElementaryOperator, ScalarOperator
-from .manipulation import OperatorArithmetics
+from typing import Any, Mapping, List, Union
+from ..operators import ElementaryOperator, OperatorArithmetics, ScalarOperator
 from .schedule import Schedule
 import warnings
 
@@ -197,16 +195,17 @@ def evaluate(
         self, op: ElementaryOperator | ScalarOperator
     ) -> CudmOperatorTerm | CallbackCoefficient | Number:
         logger.info(f"Evaluating {op}")
-        if isinstance(op, ScalarOperator):
-            if op._constant_value is None:
+        if not hasattr(op, "degrees"):
+            if not op.is_constant():
                 return CPUCallback(
-                    self._wrap_callback(op.generator, op.parameters))
+                    self._wrap_callback(lambda **kwargs: op.evaluate(**kwargs),
+                                        op.parameters))
             else:
-                return op._constant_value
+                return op.evaluate()
         else:
-            if op._id == "identity":
+            if op == op.__class__(op.degrees[0]):
                 return 1.0
-            if len(op.parameters) > 0:
+            if hasattr(op, "parameters") and len(op.parameters) > 0:
                 for param in op.parameters:
                     if param not in self._schedule._parameters:
                         raise ValueError(

python/cudaq/operator/cudm_solver.py → python/cudaq/dynamics/cudm_solver.py

@@ -11,8 +11,8 @@
 
 from .cudm_helpers import CuDensityMatOpConversion, constructLiouvillian
 from .schedule import Schedule
-from .expressions import Operator
-from cudaq.mlir._mlir_libs._quakeDialects import cudaq_runtime
+from ..operators import Operator
+from ..mlir._mlir_libs._quakeDialects import cudaq_runtime
 from .cudm_helpers import cudm, CudmStateType
 from .cudm_state import CuDensityMatState, as_cudm_state
 from .helpers import InitialState, InitialStateArgT
@@ -79,14 +79,14 @@ def evolve_dynamics(
         # Always solve the master equation if the input is a density matrix
         me_solve = True
 
-    with ScopeTimer("evolve.hamiltonian._evaluate") as timer:
-        ham_term = hamiltonian._evaluate(
+    with ScopeTimer("evolve.hamiltonian._transform") as timer:
+        ham_term = hamiltonian._transform(
             CuDensityMatOpConversion(dimensions, schedule))
     linblad_terms = []
     for c_op in collapse_operators:
-        with ScopeTimer("evolve.collapse_operators._evaluate") as timer:
+        with ScopeTimer("evolve.collapse_operators._transform") as timer:
             linblad_terms.append(
-                c_op._evaluate(CuDensityMatOpConversion(dimensions, schedule)))
+                c_op._transform(CuDensityMatOpConversion(dimensions, schedule)))
 
     with ScopeTimer("evolve.constructLiouvillian") as timer:
         liouvillian = constructLiouvillian(hilbert_space_dims, ham_term,
@@ -103,7 +103,7 @@ def evolve_dynamics(
     expectation_op = [
         cudm.Operator(
             hilbert_space_dims,
-            (observable._evaluate(CuDensityMatOpConversion(dimensions)), 1.0))
+            (observable._transform(CuDensityMatOpConversion(dimensions)), 1.0))
         for observable in observables
     ]
     integrator.set_state(initial_state, schedule._steps[0])