pyQuil v4: DefGate "matrix" property fails if unitary has mixed types

[steve-jeffrey]

Pre-Report Checklist

• [X ] I am running the latest versions of pyQuil and the Forest SDK
• [ X] I checked to make sure that this bug has not already been reported

Issue Description

A DefGate's matrix can be obtained using the "matrix" property. This property:

@property
    def matrix(self) -> np.ndarray:
        to_py_matrix = np.vectorize(_convert_to_py_expression)
        return to_py_matrix(np.asarray(super().specification.to_matrix())) 

uses numpy.vectorize to convert the gate "specification" matrix to a matrix of Expressions. The data type is set by Numpy's vectorize._vectorize_call() which doesn't set the type correctly if the matrix contains a mix of datatypes. For example, in the "works" matrix in the example shown below, the first element (quil_sin(omega)) is non-numeric, so Numpy constructs an array of objects (otype='O') and the conversion succeeds. In contrast, in the "fails" matrix in the example shown below, the first element (0) is numeric, so Numpy constructs an array of numeric types (otype='D') and then fails when it encounters the object type in the second element (quil_sin(omega)).

This problem can be overcome by setting the otype argument to np.vectorize, as shown in the example below. In the example the matrix contains Expressions, so otype='O' is required, but if the matrix only contained numeric types then otype='O' should suffice.

How to Reproduce

To demonstrate the issue, run the code as shown below.

To show how the issue can be overcome, in the run function:

• comment out the line labeled "Fails",
• enable the two lines labeled "Works".

Code Snippet

import numpy as np
from pyquil.quilatom import Parameter, quil_sin, _convert_to_py_expression
from pyquil.quilbase import DefGate

omega = Parameter("omega")

works = np.array(
    [
        [quil_sin(omega), 0],
        [0, 0],
    ]
)
fails = np.array(
    [
        [0, quil_sin(omega)],
        [0, 0],
    ]
)


def run(unitary: np.ndarray):
    gate_def = DefGate("MY_GATE", unitary, [omega])

    # Fails
    matrix = gate_def.matrix

    # Works - emulates code in the "matrix" property, but specifies the output type
    # to_py_matrix = np.vectorize(_convert_to_py_expression, otypes=['O'])
    # matrix = to_py_matrix(np.asarray(gate_def.specification.to_matrix()))

    print(f"Matrix is:\n{matrix}")


run(works)
print("-----------------------")
run(fails)

Error Output

Traceback (most recent call last):
  File "test.py", line 36, in <module>
    run(fails)
  File "test.py", line 25, in run
    matrix = gate_def.matrix
  File "...site-packages/pyquil/quilbase.py", line 682, in matrix
    return to_py_matrix(np.asarray(super().specification.to_matrix()))  # type: ignore[no-any-return]
  File "...site-packages/numpy/lib/function_base.py", line 2328, in __call__
    return self._vectorize_call(func=func, args=vargs)
  File "...site-packages/numpy/lib/function_base.py", line 2414, in _vectorize_call
    res = asanyarray(outputs, dtype=otypes[0])
TypeError: must be real number, not Function

Environment Context

pyQuil: 4.0.3

Operating System: Fedora 38

Python Version (python -V): 3.9.9