qudit-expr 0.3.4

use std::collections::{BTreeSet, HashMap};
use std::ops::{Deref, DerefMut};

use faer::Mat;
use qudit_core::ComplexScalar;
use qudit_core::QuditSystem;
use qudit_core::Radices;
use qudit_core::UnitaryMatrix;
use serde::{Deserialize, Serialize};

use crate::{ComplexExpression, UnitarySystemExpression};
use crate::{
    GenerationShape, TensorExpression,
    expressions::JittableExpression,
    index::{IndexDirection, TensorIndex},
};

use super::NamedExpression;

#[derive(PartialEq, Eq, Debug, Clone, Serialize, Deserialize)]
pub struct UnitaryExpression {
    inner: NamedExpression,
    radices: Radices,
}

// pub trait ExpressionContainer {
//     pub fn elements(&self) -> &[ComplexExpression];
//     pub fn elements_mut(&mut self) -> &mut [ComplexExpression];

//     pub fn conjugate(&mut self) {
//         todo!()
//     }
// }

// pub trait MatrixExpression: ExpressionContainer {
//     pub fn nrows(&self) -> usize;
//     pub fn ncols(&self) -> usize;

//     pub fn dagger(&mut self) {
//         todo!()
//     }

//     pub fn transpose(&mut self) {
//         todo!()
//     }
// }

impl UnitaryExpression {
    pub fn new<T: AsRef<str>>(input: T) -> Self {
        TensorExpression::new(input).try_into().unwrap()
    }

    // TODO: Change ToRadices by implementing appropriate Froms
    pub fn identity<S: Into<String>, T: Into<Radices>>(name: S, radices: T) -> Self {
        let radices = radices.into();
        let dim = radices.dimension();
        let mut body = Vec::with_capacity(dim * dim);
        for i in 0..dim {
            for j in 0..dim {
                if i == j {
                    body.push(ComplexExpression::one());
                } else {
                    body.push(ComplexExpression::zero());
                }
            }
        }
        let inner = NamedExpression::new(name, vec![], body);

        UnitaryExpression { inner, radices }
    }

    pub fn set_radices(&mut self, new_radices: Radices) {
        assert_eq!(self.radices.dimension(), new_radices.dimension());
        self.radices = new_radices;
    }

    pub fn transpose(&mut self) {
        let dim = self.dimension();
        for i in 0..dim {
            for j in (i + 1)..dim {
                let (head, tail) = self.split_at_mut(j * dim + i);
                std::mem::swap(&mut head[i * dim + j], &mut tail[0]);
            }
        }
    }

    pub fn dagger(&mut self) {
        self.conjugate();
        self.transpose();
    }

    pub fn classically_multiplex(
        expressions: &[&UnitaryExpression],
        new_dim_radices: &[usize],
    ) -> UnitarySystemExpression {
        let new_dim = new_dim_radices.iter().product();
        assert!(
            expressions.len() == new_dim,
            "Cannot multiplex a number of expressions not equal to the length of new dimension."
        );
        assert!(!expressions.is_empty());
        for expression in expressions {
            assert_eq!(
                expression.radices(),
                expressions[0].radices(),
                "All expressions must have equal radices."
            );
        }

        // construct larger tensor
        let mut new_expressions =
            Vec::with_capacity(expressions[0].dimension() * expressions[0].dimension() * new_dim);
        let mut var_map_number = 0;
        let mut new_variables = Vec::new();
        for i in 0..new_dim {
            // TODO: double clone
            let mut renamed_expression = expressions[i].clone();
            renamed_expression.alpha_rename(Some(var_map_number));
            new_expressions.extend(renamed_expression.elements().iter().cloned());
            for _ in 0..expressions[i].variables().len() {
                new_variables.push(format!("alpha_{}", var_map_number));
                var_map_number += 1;
            }
        }

        let new_indices = new_dim_radices
            .iter()
            .map(|r| (IndexDirection::Batch, (*r)))
            .chain(
                expressions[0]
                    .radices()
                    .iter()
                    .map(|r| (IndexDirection::Output, usize::from(*r))),
            )
            .chain(
                expressions[0]
                    .radices()
                    .iter()
                    .map(|r| (IndexDirection::Input, usize::from(*r))),
            )
            .enumerate()
            .map(|(id, (dir, size))| TensorIndex::new(dir, id, size))
            .collect();

        let name = {
            // if all expressions have same name: multiplex_name
            // else: multiplex_name_name_...
            if expressions
                .iter()
                .all(|e| e.name() == expressions[0].name())
            {
                format!("Multiplexed_{}", expressions[0].name())
            } else {
                "Multiplexed_".to_string()
                    + &expressions
                        .iter()
                        .map(|e| e.name())
                        .collect::<Vec<_>>()
                        .join("_")
            }
        };

        let inner = NamedExpression::new(name, new_variables, new_expressions);
        TensorExpression::from_raw(new_indices, inner)
            .try_into()
            .unwrap()
    }

    // TODO: better API for user-facing thoughts
    pub fn classically_control(
        &self,
        positions: &[usize],
        new_dim_radices: &[usize],
    ) -> UnitarySystemExpression {
        let new_dim = new_dim_radices.iter().product();
        assert!(
            positions.len() <= new_dim,
            "Cannot place unitary in more locations than length of new dimension."
        );

        // Ensure positions are unique
        let mut sorted_positions = positions.to_vec();
        sorted_positions.sort_unstable();
        assert!(
            sorted_positions.iter().collect::<BTreeSet<_>>().len() == sorted_positions.len(),
            "Positions must be unique"
        );

        // Construct identity expression
        let mut identity = Vec::with_capacity(self.dimension() * self.dimension());
        for i in 0..self.dimension() {
            for j in 0..self.dimension() {
                if i == j {
                    identity.push(ComplexExpression::one());
                } else {
                    identity.push(ComplexExpression::zero());
                }
            }
        }

        // construct larger tensor
        let mut expressions = Vec::with_capacity(self.dimension() * self.dimension() * new_dim);
        for i in 0..new_dim {
            if positions.contains(&i) {
                expressions.extend(self.elements().iter().cloned());
            } else {
                expressions.extend(identity.iter().cloned());
            }
        }

        let new_indices = new_dim_radices
            .iter()
            .map(|r| (IndexDirection::Batch, (*r)))
            .chain(
                self.radices
                    .iter()
                    .map(|r| (IndexDirection::Output, usize::from(*r))),
            )
            .chain(
                self.radices
                    .iter()
                    .map(|r| (IndexDirection::Input, usize::from(*r))),
            )
            .enumerate()
            .map(|(id, (dir, size))| TensorIndex::new(dir, id, size))
            .collect();

        let inner = NamedExpression::new(
            format!("Stacked_{}", self.name()),
            self.variables().to_owned(),
            expressions,
        );
        TensorExpression::from_raw(new_indices, inner)
            .try_into()
            .unwrap()
    }

    pub fn embed(
        &mut self,
        sub_matrix: UnitaryExpression,
        top_left_row_idx: usize,
        top_left_col_idx: usize,
    ) {
        let nrows = self.dimension();
        let ncols = self.dimension();
        let sub_nrows = sub_matrix.dimension();
        let sub_ncols = sub_matrix.dimension();

        if top_left_row_idx + sub_nrows > nrows || top_left_col_idx + sub_ncols > ncols {
            panic!("Embedding matrix is too large");
        }

        for i in 0..sub_nrows {
            for j in 0..sub_ncols {
                // TODO: remove clone by building into_iter for sub_matrix
                let sub_expr = sub_matrix[i * sub_ncols + j].clone();
                self[(top_left_row_idx + i) * ncols + (top_left_col_idx + j)] = sub_expr;
            }
        }

        // Update variables: collect all unique variables from self and sub_matrix
        let mut new_variables: Vec<String> = self.variables().to_vec();
        for var in sub_matrix.variables().iter() {
            if !new_variables.contains(var) {
                new_variables.push(var.clone());
            }
        }
        self.set_variables(new_variables);
    }

    pub fn otimes<U: AsRef<UnitaryExpression>>(&self, other: U) -> Self {
        let other = other.as_ref();
        let mut variables = Vec::new();
        let mut var_map_self = HashMap::new();
        let mut i = 0;
        for var in self.variables().iter() {
            variables.push(format!("x{}", i));
            var_map_self.insert(var.clone(), format!("x{}", i));
            i += 1;
        }
        let mut var_map_other = HashMap::new();
        for var in other.variables().iter() {
            variables.push(format!("x{}", i));
            var_map_other.insert(var.clone(), format!("x{}", i));
            i += 1;
        }

        let lhs_nrows = self.dimension();
        let lhs_ncols = self.dimension();
        let rhs_nrows = other.dimension();
        let rhs_ncols = other.dimension();

        let mut out_body = Vec::with_capacity(lhs_nrows * lhs_ncols * rhs_ncols * rhs_nrows);

        for i in 0..lhs_nrows {
            for j in 0..rhs_nrows {
                for k in 0..lhs_ncols {
                    for l in 0..rhs_ncols {
                        out_body.push(
                            self[i * lhs_ncols + k].map_var_names(&var_map_self)
                                * other[j * rhs_ncols + l].map_var_names(&var_map_other),
                        );
                    }
                }
            }
        }

        let inner = NamedExpression::new(
            format!("{} ⊗ {}", self.name(), other.name()),
            variables,
            out_body,
        );

        UnitaryExpression {
            inner,
            radices: self.radices.concat(&other.radices),
        }
    }

    pub fn dot<U: AsRef<UnitaryExpression>>(&self, other: U) -> Self {
        let other = other.as_ref();
        let lhs_nrows = self.dimension();
        let lhs_ncols = self.dimension();
        let rhs_nrows = other.dimension();
        let rhs_ncols = other.dimension();

        if lhs_ncols != rhs_nrows {
            panic!(
                "Matrix dimensions do not match for dot product: {} != {}",
                lhs_ncols, rhs_nrows
            );
        }

        let mut variables = Vec::new();
        let mut var_map_self = HashMap::new();
        let mut i = 0;
        for var in self.variables().iter() {
            variables.push(format!("x{}", i));
            var_map_self.insert(var.clone(), format!("x{}", i));
            i += 1;
        }
        let mut var_map_other = HashMap::new();
        for var in other.variables().iter() {
            variables.push(format!("x{}", i));
            var_map_other.insert(var.clone(), format!("x{}", i));
            i += 1;
        }

        let mut out_body = Vec::with_capacity(lhs_nrows * rhs_ncols);

        for i in 0..lhs_nrows {
            for j in 0..rhs_ncols {
                let mut sum = self[i * lhs_ncols].map_var_names(&var_map_self)
                    * other[j].map_var_names(&var_map_other);
                for k in 1..lhs_ncols {
                    sum += self[i * lhs_ncols + k].map_var_names(&var_map_self)
                        * other[k * rhs_ncols + j].map_var_names(&var_map_other);
                }
                out_body.push(sum);
            }
        }

        let inner = NamedExpression::new(
            format!("{} ⋅ {}", self.name(), other.name()),
            variables,
            out_body,
        );

        UnitaryExpression {
            inner,
            radices: self.radices.clone(),
        }
    }

    pub fn get_arg_map<C: ComplexScalar>(&self, args: &[C::R]) -> HashMap<&str, C::R> {
        self.variables()
            .iter()
            .zip(args.iter())
            .map(|(a, b)| (a.as_str(), *b))
            .collect()
    }

    pub fn eval<C: ComplexScalar>(&self, args: &[C::R]) -> UnitaryMatrix<C> {
        let arg_map = self.get_arg_map::<C>(args);
        let dim = self.radices.dimension();
        let mut mat = Mat::zeros(dim, dim);
        for i in 0..dim {
            for j in 0..dim {
                *mat.get_mut(i, j) = self[i * self.dimension() + j].eval(&arg_map);
            }
        }
        UnitaryMatrix::new(self.radices.clone(), mat)
    }

    /// Evaluates the matrix-by-vector derivative of this expression at `args`.
    ///
    /// Returns one matrix per parameter (in the same order as `variables()`),
    /// where the `k`-th matrix is the elementwise partial derivative of the
    /// unitary with respect to the `k`-th parameter.
    pub fn eval_grad<C: ComplexScalar>(&self, args: &[C::R]) -> Vec<Mat<C>> {
        let arg_map = self.get_arg_map::<C>(args);
        let dim = self.radices.dimension();
        self.variables()
            .iter()
            .map(|var| {
                let mut mat = Mat::zeros(dim, dim);
                for i in 0..dim {
                    for j in 0..dim {
                        *mat.get_mut(i, j) = self[i * dim + j].differentiate(var).eval(&arg_map);
                    }
                }
                mat
            })
            .collect()
    }

    /// Returns the canonical OpenQASM 2.0 (`qelib1.inc`) gate name for this
    /// expression, if one exists, based on its name, qudit count, and
    /// parameter count. Returns `None` for expressions with no fixed QASM
    /// 2.0 equivalent (e.g. arbitrary/custom unitaries).
    pub fn qasm_name(&self) -> Option<&'static str> {
        qasm_gate_name_lookup(self.name(), self.radices().len(), self.num_params())
    }
}

pub(crate) fn qasm_gate_name_lookup(
    name: &str,
    n_qudits: usize,
    n_params: usize,
) -> Option<&'static str> {
    let mut name = name;
    while let Some(stripped) = name.strip_suffix("_subbed") {
        name = stripped;
    }
    name = name.strip_prefix("Stacked_").unwrap_or(name);

    match (name, n_qudits, n_params) {
        // Built-in primitives
        ("U", 1, 3) => Some("U"),
        // qelib1.inc standard gates
        ("U2", 1, 2) => Some("u2"),
        ("U1", 1, 1) => Some("u1"),
        ("I", 1, 0) => Some("id"),
        ("X", 1, 0) => Some("x"),
        ("Y", 1, 0) => Some("y"),
        ("Z", 1, 0) => Some("z"),
        ("H", 1, 0) => Some("h"),
        ("S", 1, 0) => Some("s"),
        ("T", 1, 0) => Some("t"),
        ("SX", 1, 0) => Some("sx"),
        ("P", 1, 1) => Some("p"),
        ("RX", 1, 1) => Some("rx"),
        ("RY", 1, 1) => Some("ry"),
        ("RZ", 1, 1) => Some("rz"),
        ("Swap", 2, 0) => Some("swap"),
        ("RXX", 2, 1) => Some("rxx"),
        ("RZZ", 2, 1) => Some("rzz"),
        ("Controlled(X)", 2, 0) => Some("cx"),
        ("Controlled(X)", 3, 0) => Some("ccx"),
        ("Controlled(X)", 4, 0) => Some("c3x"),
        ("Controlled(X)", 5, 0) => Some("c4x"),
        ("Controlled(Z)", 2, 0) => Some("cz"),
        ("Controlled(Y)", 2, 0) => Some("cy"),
        ("Controlled(H)", 2, 0) => Some("ch"),
        ("Controlled(Swap)", 3, 0) => Some("cswap"),
        ("Controlled(SX)", 2, 0) => Some("csx"),
        ("Controlled(RX)", 2, 1) => Some("crx"),
        ("Controlled(RY)", 2, 1) => Some("cry"),
        ("Controlled(RZ)", 2, 1) => Some("crz"),
        ("Controlled(U1)", 2, 1) => Some("cu1"),
        ("Controlled(P)", 2, 1) => Some("cp"),
        ("Controlled(U)", 2, 3) => Some("cu3"),
        ("Dagger(S)", 1, 0) => Some("sdg"),
        ("Dagger(T)", 1, 0) => Some("tdg"),
        ("Dagger(SX)", 1, 0) => Some("sxdg"),
        _ => None,
    }
}

impl JittableExpression for UnitaryExpression {
    fn generation_shape(&self) -> GenerationShape {
        GenerationShape::Matrix(self.radices.dimension(), self.radices.dimension())
    }
}

impl AsRef<UnitaryExpression> for UnitaryExpression {
    fn as_ref(&self) -> &UnitaryExpression {
        self
    }
}

impl AsRef<NamedExpression> for UnitaryExpression {
    fn as_ref(&self) -> &NamedExpression {
        &self.inner
    }
}

impl From<UnitaryExpression> for NamedExpression {
    fn from(value: UnitaryExpression) -> Self {
        value.inner
    }
}

impl Deref for UnitaryExpression {
    type Target = NamedExpression;

    fn deref(&self) -> &Self::Target {
        &self.inner
    }
}

impl DerefMut for UnitaryExpression {
    fn deref_mut(&mut self) -> &mut Self::Target {
        &mut self.inner
    }
}

impl From<UnitaryExpression> for TensorExpression {
    fn from(value: UnitaryExpression) -> Self {
        let UnitaryExpression { inner, radices } = value;
        let indices = radices
            .iter()
            .map(|r| (IndexDirection::Output, usize::from(*r)))
            .chain(
                radices
                    .iter()
                    .map(|r| (IndexDirection::Input, usize::from(*r))),
            )
            .enumerate()
            .map(|(i, (d, r))| TensorIndex::new(d, i, r))
            .collect();
        TensorExpression::from_raw(indices, inner)
    }
}

impl TryFrom<TensorExpression> for UnitaryExpression {
    // TODO: Come up with proper error handling
    type Error = String;

    fn try_from(value: TensorExpression) -> Result<Self, Self::Error> {
        let mut input_radices = vec![];
        let mut output_radices = vec![];
        for idx in value.indices() {
            match idx.direction() {
                IndexDirection::Input => {
                    input_radices.push(idx.index_size());
                }
                IndexDirection::Output => {
                    output_radices.push(idx.index_size());
                }
                _ => {
                    return Err(String::from(
                        "Cannot convert a tensor with non-input, non-output indices to an isometry.",
                    ));
                }
            }
        }

        if input_radices != output_radices {
            return Err(String::from(
                "Non-square matrix tensor cannot be converted to a unitary.",
            ));
        }

        Ok(UnitaryExpression {
            inner: value.into(),
            radices: input_radices.into(),
        })
    }
}

impl<C: ComplexScalar> From<UnitaryMatrix<C>> for UnitaryExpression {
    fn from(value: UnitaryMatrix<C>) -> Self {
        let mut body = vec![Vec::with_capacity(value.ncols()); value.nrows()];
        for col in value.col_iter() {
            for (row_id, elem) in col.iter().enumerate() {
                body[row_id].push(ComplexExpression::from(*elem));
            }
        }
        let mut flat_body = Vec::with_capacity(value.ncols() * value.nrows());
        for row in body.into_iter() {
            for elem in row.into_iter() {
                flat_body.push(elem);
            }
        }
        let inner = NamedExpression::new("Constant", vec![], flat_body);

        UnitaryExpression {
            inner,
            radices: value.radices(),
        }
    }
}

impl QuditSystem for UnitaryExpression {
    fn radices(&self) -> Radices {
        self.radices.clone()
    }
}

#[cfg(feature = "python")]
mod python {
    use super::*;
    use crate::python::PyExpressionRegistrar;
    use numpy::PyArray2;
    use numpy::PyArray3;
    use numpy::PyArrayMethods;
    use numpy::ndarray::{ArrayViewMut2, ArrayViewMut3};
    use pyo3::{
        PyResult,
        exceptions::PyTypeError,
        prelude::*,
        types::{PyBytes, PyTuple},
    };
    use pyo3_stub_gen::derive::*;
    use pyo3_stub_gen::impl_stub_type;
    use qudit_core::Radix;
    use qudit_core::c64;

    /// A symbolic, parameterized expression representing a unitary matrix.
    #[gen_stub_pyclass]
    #[pyclass(name = "UnitaryExpression", module = "openqudit.expressions")]
    pub struct PyUnitaryExpression {
        expr: UnitaryExpression,
    }

    #[gen_stub_pymethods]
    #[pymethods]
    impl PyUnitaryExpression {
        /// Parses a unitary expression from its string representation.
        ///
        /// # Arguments
        ///
        /// * `expr` - The textual definition of the unitary expression.
        #[new]
        fn new(expr: String) -> Self {
            Self {
                expr: UnitaryExpression::new(expr),
            }
        }

        /// Constructs an identity unitary expression over the given radices.
        ///
        /// # Arguments
        ///
        /// * `name` - The name to assign to the resulting expression.
        /// * `radices` - The radix of each qudit in the system.
        #[staticmethod]
        fn identity(name: String, radices: Vec<usize>) -> Self {
            Self {
                expr: UnitaryExpression::identity(name, radices),
            }
        }

        /// Evaluates this expression at the given parameter values and returns
        /// the resulting unitary matrix as a NumPy array.
        ///
        /// # Arguments
        ///
        /// * `args` - The real-valued parameters to substitute into the expression,
        ///   in the same order as `variables()`.
        #[pyo3(signature = (*args))]
        fn __call__<'py>(&self, args: &Bound<'py, PyTuple>) -> PyResult<Bound<'py, PyArray2<c64>>> {
            let py = args.py();
            let args: Vec<f64> = args.extract()?;
            let unitary = self.expr.eval(&args);
            let py_array: Bound<'py, PyArray2<c64>> =
                PyArray2::zeros(py, (unitary.dimension(), unitary.dimension()), false);

            {
                let mut readwrite = py_array.readwrite();
                let mut py_array_view: ArrayViewMut2<c64> = readwrite.as_array_mut();

                for (j, col) in unitary.col_iter().enumerate() {
                    for (i, val) in col.iter().enumerate() {
                        py_array_view[[i, j]] = *val;
                    }
                }
            }

            Ok(py_array)
        }

        /// Evaluates the matrix-by-vector derivative of this expression at the
        /// given parameter values.
        ///
        /// Returns a NumPy array of shape `(num_params, dim, dim)`, where the
        /// `k`-th entry along the first axis is the elementwise partial
        /// derivative of the unitary matrix with respect to the `k`-th
        /// parameter, in the same order as `variables()`.
        ///
        /// # Arguments
        ///
        /// * `args` - The real-valued parameters to substitute into the expression,
        ///   in the same order as `variables()`.
        #[pyo3(signature = (*args))]
        fn gradient<'py>(&self, args: &Bound<'py, PyTuple>) -> PyResult<Bound<'py, PyArray3<c64>>> {
            let py = args.py();
            let args: Vec<f64> = args.extract()?;
            let dim = self.expr.dimension();
            let num_params = self.expr.num_params();
            let grad = self.expr.eval_grad(&args);
            let py_array: Bound<'py, PyArray3<c64>> =
                PyArray3::zeros(py, (num_params, dim, dim), false);

            {
                let mut readwrite = py_array.readwrite();
                let mut py_array_view: ArrayViewMut3<c64> = readwrite.as_array_mut();

                for (p, mat) in grad.iter().enumerate() {
                    for (j, col) in mat.col_iter().enumerate() {
                        for (i, val) in col.iter().enumerate() {
                            py_array_view[[p, i, j]] = *val;
                        }
                    }
                }
            }

            Ok(py_array)
        }

        /// Returns the number of free (unbound) parameters in this expression.
        fn num_params(&self) -> usize {
            self.expr.num_params()
        }

        /// Returns the name assigned to this expression.
        fn name(&self) -> String {
            self.expr.name().to_string()
        }

        /// Returns the radix of each qudit that this unitary acts on.
        fn radices(&self) -> Vec<Radix> {
            self.expr.radices().to_vec()
        }

        /// Returns the total Hilbert space dimension of the underlying qudit system.
        fn dimension(&self) -> usize {
            self.expr.dimension()
        }

        /// Returns the canonical OpenQASM 2.0 (`qelib1.inc`) gate name for this
        /// expression (e.g. `"rz"`, `"cx"`), or `None` if it has no fixed QASM
        /// 2.0 equivalent, such as an arbitrary or custom unitary.
        fn qasm_name(&self) -> Option<&'static str> {
            self.expr.qasm_name()
        }

        /// Transposes this unitary expression in place.
        fn transpose(&mut self) {
            self.expr.transpose();
        }

        /// Conjugate-transposes (Hermitian adjoint) this unitary expression in place.
        fn dagger(&mut self) {
            self.expr.dagger();
        }

        /// Computes the tensor (Kronecker) product of this expression with `other`,
        /// returning a new expression over the combined qudit system.
        ///
        /// # Arguments
        ///
        /// * `other` - The unitary expression to tensor with this one.
        fn otimes(&self, other: &PyUnitaryExpression) -> Self {
            Self {
                expr: self.expr.otimes(&other.expr),
            }
        }

        /// Computes the matrix product of this expression with `other`,
        /// returning a new expression.
        ///
        /// # Arguments
        ///
        /// * `other` - The unitary expression to multiply with this one.
        fn dot(&self, other: &PyUnitaryExpression) -> Self {
            Self {
                expr: self.expr.dot(&other.expr),
            }
        }

        /// Embeds `sub_matrix` into this expression's matrix in place, placing its
        /// top-left corner at the given row and column index.
        ///
        /// # Arguments
        ///
        /// * `sub_matrix` - The smaller unitary expression to embed.
        /// * `top_left_row_idx` - Row index at which to place the sub-matrix.
        /// * `top_left_col_idx` - Column index at which to place the sub-matrix.
        fn embed(
            &mut self,
            sub_matrix: &PyUnitaryExpression,
            top_left_row_idx: usize,
            top_left_col_idx: usize,
        ) {
            self.expr
                .embed(sub_matrix.expr.clone(), top_left_row_idx, top_left_col_idx);
        }

        // fn classically_control(&self, positions: Vec<usize>, new_dim_radices: Vec<usize>) -> crate::python::tensor::PyTensorExpression {
        //     let result = self.expr.classically_control(&positions, &new_dim_radices);
        //     result.into()
        // }

        fn __repr__(&self) -> String {
            format!(
                "UnitaryExpression(name='{}', radices={:?}, params={})",
                self.expr.name(),
                self.expr.radices().to_vec(),
                self.expr.num_params()
            )
        }

        fn __matmul__(&self, other: UnitaryExpression) -> PyUnitaryExpression {
            Self {
                expr: self.expr.dot(&other),
            }
        }

        pub fn __setstate__(&mut self, state: &Bound<'_, PyBytes>) -> PyResult<()> {
            self.expr = postcard::from_bytes(state.as_bytes())
                .map_err(|e| PyTypeError::new_err(format!("Failed to deserialize circuit: {e}")))?;
            Ok(())
        }

        pub fn __getstate__<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyBytes>> {
            let bytes: Vec<u8> = postcard::to_allocvec(&self.expr)
                .map_err(|e| PyTypeError::new_err(format!("Failed to serialize circuit: {e}")))?;
            Ok(PyBytes::new(py, &bytes))
        }

        pub fn __getnewargs__(&self) -> PyResult<(String,)> {
            // Dummy arg: __setstate__ fully restores the expression, so we just
            // need any valid value for the expr argument of __new__.
            Ok(("I<2>(){[[1,0,],[0,1,],]}".to_string(),))
        }
    }

    impl From<UnitaryExpression> for PyUnitaryExpression {
        fn from(value: UnitaryExpression) -> Self {
            PyUnitaryExpression { expr: value }
        }
    }

    impl From<PyUnitaryExpression> for UnitaryExpression {
        fn from(value: PyUnitaryExpression) -> Self {
            value.expr
        }
    }

    impl_stub_type!(UnitaryExpression = PyUnitaryExpression);

    impl<'py> IntoPyObject<'py> for UnitaryExpression {
        type Target = <PyUnitaryExpression as IntoPyObject<'py>>::Target;
        type Output = Bound<'py, Self::Target>;
        type Error = PyErr;

        fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
            let py_expr = PyUnitaryExpression::from(self);
            Bound::new(py, py_expr)
        }
    }

    impl<'a, 'py> FromPyObject<'a, 'py> for UnitaryExpression {
        type Error = PyErr;

        fn extract(ob: Borrowed<'a, 'py, PyAny>) -> PyResult<Self> {
            let py_expr: PyRef<PyUnitaryExpression> = ob.extract()?;
            Ok(py_expr.expr.clone())
        }
    }

    /// Registers the UnitaryExpression class with the Python module.
    fn register(parent_module: &Bound<'_, PyModule>) -> PyResult<()> {
        parent_module.add_class::<PyUnitaryExpression>()?;
        Ok(())
    }
    inventory::submit!(PyExpressionRegistrar { func: register });
}