diff --git a/qualtran/_infra/data_types.py b/qualtran/_infra/data_types.py index 227247191..015a3901f 100644 --- a/qualtran/_infra/data_types.py +++ b/qualtran/_infra/data_types.py @@ -77,10 +77,30 @@ def get_classical_domain(self) -> Iterable[Any]: def to_bits(self, x) -> List[int]: """Yields individual bits corresponding to binary representation of x""" + def to_bits_array(self, x_array: NDArray[Any]) -> NDArray[np.uint8]: + """Yields an NDArray of bits corresponding to binary representations of the input elements. + + Often, converting an array can be performed faster than converting each element individually. + This operation accepts any NDArray of values, and the output array satisfies + `output_shape = input_shape + (self.bitsize,)`. + """ + return np.vectorize( + lambda x: np.asarray(self.to_bits(x), dtype=np.uint8), signature='()->(n)' + )(x_array) + @abc.abstractmethod def from_bits(self, bits: Sequence[int]): """Combine individual bits to form x""" + def from_bits_array(self, bits_array: NDArray[np.uint8]): + """Combine individual bits to form classical values. + + Often, converting an array can be performed faster than converting each element individually. + This operation accepts any NDArray of bits such that the last dimension equals `self.bitsize`, + and the output array satisfies `output_shape = input_shape[:-1]`. + """ + return np.vectorize(self.from_bits, signature='(n)->()')(bits_array) + @abc.abstractmethod def assert_valid_classical_val(self, val: Any, debug_str: str = 'val'): """Raises an exception if `val` is not a valid classical value for this type. @@ -90,17 +110,6 @@ def assert_valid_classical_val(self, val: Any, debug_str: str = 'val'): debug_str: Optional debugging information to use in exception messages. """ - @abc.abstractmethod - def is_symbolic(self) -> bool: - """Returns True if this qdtype is parameterized with symbolic objects.""" - - def iteration_length_or_zero(self) -> SymbolicInt: - """Safe version of iteration length. - - Returns the iteration_length if the type has it or else zero. - """ - return getattr(self, 'iteration_length', 0) - def assert_valid_classical_val_array(self, val_array: NDArray[Any], debug_str: str = 'val'): """Raises an exception if `val_array` is not a valid array of classical values for this type. @@ -116,6 +125,17 @@ def assert_valid_classical_val_array(self, val_array: NDArray[Any], debug_str: s for val in val_array.reshape(-1): self.assert_valid_classical_val(val) + @abc.abstractmethod + def is_symbolic(self) -> bool: + """Returns True if this qdtype is parameterized with symbolic objects.""" + + def iteration_length_or_zero(self) -> SymbolicInt: + """Safe version of iteration length. + + Returns the iteration_length if the type has it or else zero. + """ + return getattr(self, 'iteration_length', 0) + def __str__(self): return f'{self.__class__.__name__}({self.num_qubits})' @@ -324,10 +344,43 @@ def to_bits(self, x: int) -> List[int]: self.assert_valid_classical_val(x) return [int(x) for x in f'{int(x):0{self.bitsize}b}'] + def to_bits_array(self, x_array: NDArray[np.integer]) -> NDArray[np.uint8]: + """Returns the big-endian bitstrings specified by the given integers. + + Args: + x_array: An integer or array of unsigned integers. + """ + if is_symbolic(self.bitsize): + raise ValueError(f"Cannot compute bits for symbolic {self.bitsize=}") + + w = int(self.bitsize) + x = np.atleast_1d(x_array) + if not np.issubdtype(x.dtype, np.uint): + assert np.all(x >= 0) + assert np.iinfo(x.dtype).bits <= 64 + x = x.astype(np.uint64) + assert w <= np.iinfo(x.dtype).bits + mask = 2 ** np.arange(w - 1, 0 - 1, -1, dtype=x.dtype).reshape((w, 1)) + return (x & mask).astype(bool).astype(np.uint8).T + def from_bits(self, bits: Sequence[int]) -> int: """Combine individual bits to form x""" return int("".join(str(x) for x in bits), 2) + def from_bits_array(self, bits_array: NDArray[np.uint8]) -> NDArray[np.integer]: + """Returns the integer specified by the given big-endian bitstrings. + + Args: + bits_array: A bitstring or array of bitstrings, each of which has the 1s bit (LSB) at the end. + Returns: + An array of integers; one for each bitstring. + """ + bitstrings = np.atleast_2d(bits_array) + if bitstrings.shape[1] > 64: + raise NotImplementedError() + basis = 2 ** np.arange(bitstrings.shape[1] - 1, 0 - 1, -1, dtype=np.uint64) + return np.sum(basis * bitstrings, axis=1, dtype=np.uint64) + def assert_valid_classical_val(self, val: int, debug_str: str = 'val'): if not isinstance(val, (int, np.integer)): raise ValueError(f"{debug_str} should be an integer, not {val!r}") @@ -486,12 +539,31 @@ def num_int(self) -> SymbolicInt: return self.bitsize - self.num_frac - int(self.signed) @property - def fxp_dtype_str(self) -> str: - return f'fxp-{"us"[self.signed]}{self.bitsize}/{self.num_frac}' + def fxp_dtype_template(self) -> Fxp: + """A template of the `Fxp` data type for classical values. + + - op_sizing='same' and const_op_sizing='same' ensure that the returned object is not resized + to a bigger fixed point number when doing operations with other Fxp objects. + - shifting='trunc' ensures that when shifting the Fxp integer to left / right; the digits are + truncated and no rounding occurs + - overflow='wrap' ensures that when performing operations where result overflows, the overflowed + digits are simply discarded. + """ + if is_symbolic(self.bitsize) or is_symbolic(self.num_frac): + raise ValueError( + "Cannot construct Fxp template for symbolic bitsizes {self.bitsize=}, {self.num_frac=}" + ) - @property - def _fxp_dtype(self) -> Fxp: - return Fxp(None, dtype=self.fxp_dtype_str) + return Fxp( + None, + n_word=self.bitsize, + n_frac=self.num_frac, + signed=self.signed, + op_sizing='same', + const_op_sizing='same', + shifting='trunc', + overflow='wrap', + ) def is_symbolic(self) -> bool: return is_symbolic(self.bitsize, self.num_frac) @@ -517,7 +589,7 @@ def to_bits( sign = int(x < 0) x = abs(x) fxp = x if isinstance(x, Fxp) else Fxp(x) - bits = [int(x) for x in fxp.like(self._fxp_dtype).bin()] + bits = [int(x) for x in fxp.like(self.fxp_dtype_template).bin()] if self.signed and not complement: bits[0] = sign return bits @@ -526,7 +598,14 @@ def from_bits(self, bits: Sequence[int]) -> Fxp: """Combine individual bits to form x""" bits_bin = "".join(str(x) for x in bits[:]) fxp_bin = "0b" + bits_bin[: -self.num_frac] + "." + bits_bin[-self.num_frac :] - return Fxp(fxp_bin, dtype=self.fxp_dtype_str) + return Fxp(fxp_bin, like=self.fxp_dtype_template) + + def from_bits_array(self, bits_array: NDArray[np.uint8]): + assert isinstance(self.bitsize, int), "cannot convert to bits for symbolic bitsize" + # TODO figure out why `np.vectorize` is not working here + return Fxp( + [self.from_bits(bitstring) for bitstring in bits_array.reshape(-1, self.bitsize)] + ) def to_fixed_width_int(self, x: Union[float, Fxp]) -> int: """Returns the interpretation of the binary representation of `x` as an integer. Requires `x` to be nonnegative.""" @@ -546,19 +625,37 @@ def __attrs_post_init__(self): def get_classical_domain(self) -> Iterable[Fxp]: qint = QIntOnesComp(self.bitsize) if self.signed else QUInt(self.bitsize) for x in qint.get_classical_domain(): - yield Fxp(x / 2**self.num_frac, dtype=self.fxp_dtype_str) + yield Fxp(x / 2**self.num_frac).like(self.fxp_dtype_template) def _assert_valid_classical_val(self, val: Union[float, Fxp], debug_str: str = 'val'): fxp_val = val if isinstance(val, Fxp) else Fxp(val) - if fxp_val.get_val() != fxp_val.like(self._fxp_dtype).get_val(): + if fxp_val.get_val() != fxp_val.like(self.fxp_dtype_template).get_val(): raise ValueError( f"{debug_str}={val} cannot be accurately represented using Fxp {fxp_val}" ) def assert_valid_classical_val(self, val: Union[float, Fxp], debug_str: str = 'val'): - # TODO: Asserting a valid value here opens a can of worms because classical data, except integers, - # is currently not propagated correctly through Bloqs - pass + assert isinstance(val, Fxp) + assert val.overflow == 'wrap' + assert val.shifting == 'trunc' + self._assert_valid_classical_val(val, debug_str) + + def float_to_fxp( + self, val: Union[float, int], *, raw: bool = False, require_exact: bool = True + ) -> Fxp: + """Convert a floating point value to an Fxp constant of this dtype. + + If `raw` is True, then returns `val / 2**self.n_frac` instead. + + Args: + val: Floating point value. + raw: Convert from a raw integer value instead + require_exact: If True, represent the input `val` exactly and raise + a ValueError if it cannot be represented. + """ + if require_exact: + self._assert_valid_classical_val(val if not raw else val / 2**self.num_frac) + return Fxp(val, raw=raw, like=self.fxp_dtype_template) def __str__(self): if self.signed: diff --git a/qualtran/_infra/data_types_test.py b/qualtran/_infra/data_types_test.py index dbe063045..20936eef3 100644 --- a/qualtran/_infra/data_types_test.py +++ b/qualtran/_infra/data_types_test.py @@ -11,13 +11,15 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. - import math import random +from typing import Any, Sequence, Union +import cirq import numpy as np import pytest import sympy +from numpy.typing import NDArray from qualtran.symbolics import is_symbolic @@ -95,12 +97,12 @@ def test_qfxp(): assert str(qfp_16) == 'QFxp(16, 15)' assert qfp_16.num_qubits == 16 assert qfp_16.num_int == 1 - assert qfp_16.fxp_dtype_str == 'fxp-u16/15' + assert qfp_16.fxp_dtype_template.dtype == 'fxp-u16/15' qfp_16 = QFxp(16, 15, signed=True) assert str(qfp_16) == 'QFxp(16, 15, True)' assert qfp_16.num_qubits == 16 assert qfp_16.num_int == 0 - assert qfp_16.fxp_dtype_str == 'fxp-s16/15' + assert qfp_16.fxp_dtype_template.dtype == 'fxp-s16/15' with pytest.raises(ValueError, match="num_qubits must be > 1."): QFxp(1, 1, signed=True) QFxp(1, 1, signed=False) @@ -233,8 +235,20 @@ def test_single_qubit_consistency(): assert check_dtypes_consistent(QFxp(1, 1), QBit()) -def test_to_and_from_bits(): - # QInt +def assert_to_and_from_bits_array_consistent(qdtype: QDType, values: Union[Sequence[Any], NDArray]): + values = np.asarray(values) + bits_array = qdtype.to_bits_array(values) + + # individual values + for val, bits in zip(values.reshape(-1), bits_array.reshape(-1, qdtype.num_qubits)): + assert np.all(bits == qdtype.to_bits(val)) + + # round trip + values_roundtrip = qdtype.from_bits_array(bits_array) + assert np.all(values_roundtrip == values) + + +def test_qint_to_and_from_bits(): qint4 = QInt(4) assert [*qint4.get_classical_domain()] == [*range(-8, 8)] for x in range(-8, 8): @@ -246,7 +260,10 @@ def test_to_and_from_bits(): with pytest.raises(ValueError): QInt(4).to_bits(10) - # QUInt + assert_to_and_from_bits_array_consistent(qint4, range(-8, 8)) + + +def test_quint_to_and_from_bits(): quint4 = QUInt(4) assert [*quint4.get_classical_domain()] == [*range(0, 16)] assert list(quint4.to_bits(10)) == [1, 0, 1, 0] @@ -259,23 +276,66 @@ def test_to_and_from_bits(): with pytest.raises(ValueError): quint4.to_bits(-1) - # BoundedQUInt + assert_to_and_from_bits_array_consistent(quint4, range(0, 16)) + + +def test_bits_to_int(): + rs = np.random.RandomState(52) + bitstrings = rs.choice([0, 1], size=(100, 23)) + + nums = QUInt(23).from_bits_array(bitstrings) + assert nums.shape == (100,) + + for num, bs in zip(nums, bitstrings): + ref_num = cirq.big_endian_bits_to_int(bs.tolist()) + assert num == ref_num + + # check one input bitstring instead of array of input bitstrings. + (num,) = QUInt(23).from_bits_array(np.array([1, 0])) + assert num == 2 + + +def test_int_to_bits(): + rs = np.random.RandomState(52) + nums = rs.randint(0, 2**23 - 1, size=(100,), dtype=np.uint64) + bitstrings = QUInt(23).to_bits_array(nums) + assert bitstrings.shape == (100, 23) + + for num, bs in zip(nums, bitstrings): + ref_bs = cirq.big_endian_int_to_bits(int(num), bit_count=23) + np.testing.assert_array_equal(ref_bs, bs) + + # check bounds + with pytest.raises(AssertionError): + QUInt(8).to_bits_array(np.array([4, -2])) + + +def test_bounded_quint_to_and_from_bits(): bquint4 = BoundedQUInt(4, 12) assert [*bquint4.get_classical_domain()] == [*range(0, 12)] assert list(bquint4.to_bits(10)) == [1, 0, 1, 0] with pytest.raises(ValueError): BoundedQUInt(4, 12).to_bits(13) - # QBit + assert_to_and_from_bits_array_consistent(bquint4, range(0, 12)) + + +def test_qbit_to_and_from_bits(): assert list(QBit().to_bits(0)) == [0] assert list(QBit().to_bits(1)) == [1] with pytest.raises(ValueError): QBit().to_bits(2) - # QAny + assert_to_and_from_bits_array_consistent(QBit(), [0, 1]) + + +def test_qany_to_and_from_bits(): assert list(QAny(4).to_bits(10)) == [1, 0, 1, 0] - # QIntOnesComp + assert_to_and_from_bits_array_consistent(QAny(4), range(16)) + + +def test_qintonescomp_to_and_from_bits(): qintones4 = QIntOnesComp(4) assert list(qintones4.to_bits(-2)) == [1, 1, 0, 1] assert list(qintones4.to_bits(2)) == [0, 0, 1, 0] @@ -287,6 +347,10 @@ def test_to_and_from_bits(): with pytest.raises(ValueError): qintones4.to_bits(-8) + assert_to_and_from_bits_array_consistent(qintones4, range(-7, 8)) + + +def test_qfxp_to_and_from_bits(): # QFxp: Negative numbers are stored as ones complement qfxp_4_3 = QFxp(4, 3, True) assert list(qfxp_4_3.to_bits(0.5)) == [0, 1, 0, 0] @@ -321,6 +385,11 @@ def test_to_and_from_bits(): assert list(QFxp(7, 3, True).to_bits(-4.375)) == [1] + [0, 1, 1] + [1, 0, 1] assert list(QFxp(7, 3, True).to_bits(+4.625)) == [0] + [1, 0, 0] + [1, 0, 1] + assert_to_and_from_bits_array_consistent(QFxp(4, 3, False), [1 / 2, 1 / 4, 3 / 8]) + assert_to_and_from_bits_array_consistent( + QFxp(4, 3, True), [1 / 2, -1 / 2, 1 / 4, -1 / 4, -3 / 8, 3 / 8] + ) + def test_iter_bits(): assert QUInt(2).to_bits(0) == [0, 0] diff --git a/qualtran/bloqs/arithmetic/multiplication.py b/qualtran/bloqs/arithmetic/multiplication.py index 33b85e712..8d8eeaa98 100644 --- a/qualtran/bloqs/arithmetic/multiplication.py +++ b/qualtran/bloqs/arithmetic/multiplication.py @@ -17,6 +17,7 @@ import cirq import numpy as np from attrs import evolve, frozen +from fxpmath import Fxp from qualtran import ( Bloq, @@ -57,11 +58,13 @@ def pretty_name(self) -> str: @property def signature(self) -> 'Signature': return Signature.build_from_dtypes( - a=QUInt(self.a_bitsize), - b=QUInt(self.b_bitsize), - result=QFxp(self.result_bitsize, self.result_bitsize), + a=QUInt(self.a_bitsize), b=QUInt(self.b_bitsize), result=self.result_dtype ) + @property + def result_dtype(self): + return QFxp(self.result_bitsize, self.result_bitsize) + def registers(self) -> Sequence[Union[int, Sequence[int]]]: if not isinstance(self.a_bitsize, int): raise ValueError(f'Symbolic bitsize {self.a_bitsize} not supported') @@ -80,8 +83,11 @@ def apply(self, a: int, b: int, result: int) -> Union[int, Iterable[int]]: def with_registers(self, *new_registers: Union[int, Sequence[int]]): raise NotImplementedError("Not needed.") - def on_classical_vals(self, a: int, b: int, result: int) -> Dict[str, 'ClassicalValT']: - result_out = (result + a * b * ((-1) ** self.is_adjoint)) % (2**self.result_bitsize) + def on_classical_vals(self, a: int, b: int, result: Fxp) -> Dict[str, 'ClassicalValT']: + result_out_raw = (int(result.raw()) + a * b * ((-1) ** self.is_adjoint)) % ( + 2**self.result_bitsize + ) + result_out = self.result_dtype.float_to_fxp(result_out_raw, raw=True) return {'a': a, 'b': b, 'result': result_out} def _circuit_diagram_info_(self, args: cirq.CircuitDiagramInfoArgs) -> cirq.CircuitDiagramInfo: diff --git a/qualtran/bloqs/arithmetic/multiplication_test.py b/qualtran/bloqs/arithmetic/multiplication_test.py index 791182778..2bb62ed08 100644 --- a/qualtran/bloqs/arithmetic/multiplication_test.py +++ b/qualtran/bloqs/arithmetic/multiplication_test.py @@ -159,14 +159,19 @@ def test_plus_equal_product(): for b in range(2**b_bit): for result in range(2**res_bit): res_out = (result + a * b) % 2**res_bit + # Test Bloq style classical simulation. - assert bloq.call_classically(a=a, b=b, result=result) == (a, b, res_out) + result_fxp = bloq.result_dtype.float_to_fxp(result, raw=True) + res_out_fxp = bloq.result_dtype.float_to_fxp(res_out, raw=True) + assert bloq.call_classically(a=a, b=b, result=result_fxp) == (a, b, res_out_fxp) + # Prepare basis states mapping for cirq-style simulation. input_state_str = f'{a:0{a_bit}b}' + f'{b:0{b_bit}b}' + f'{result:0{res_bit}b}' input_state = int(input_state_str, 2) output_state_str = f'{a:0{a_bit}b}' + f'{b:0{b_bit}b}' + f'{res_out:0{res_bit}b}' output_state = int(output_state_str, 2) basis_map[input_state] = output_state + # Test cirq style simulation. assert len(basis_map) == len(set(basis_map.values())) circuit = cirq.Circuit(bloq.on(*cirq.LineQubit.range(num_bits))) diff --git a/qualtran/bloqs/basic_gates/z_basis.py b/qualtran/bloqs/basic_gates/z_basis.py index bef3d03af..930a49dde 100644 --- a/qualtran/bloqs/basic_gates/z_basis.py +++ b/qualtran/bloqs/basic_gates/z_basis.py @@ -33,6 +33,7 @@ DecomposeTypeError, QAny, QBit, + QDType, Register, Side, Signature, @@ -42,7 +43,6 @@ from qualtran.bloqs.bookkeeping import ArbitraryClifford from qualtran.cirq_interop.t_complexity_protocol import TComplexity from qualtran.drawing import Circle, directional_text_box, Text, TextBox, WireSymbol -from qualtran.simulation.classical_sim import ints_to_bits if TYPE_CHECKING: import cirq @@ -384,12 +384,16 @@ def check(self, attribute, val): if val >= 2**self.bitsize: raise ValueError(f"`val` is too big for bitsize {self.bitsize}") + @cached_property + def dtype(self) -> QDType: + if self.bitsize == 1: + return QBit() + return QAny(self.bitsize) + @cached_property def signature(self) -> Signature: side = Side.RIGHT if self.state else Side.LEFT - if self.bitsize == 1: - return Signature([Register('val', QBit(), side=side)]) - return Signature([Register('val', QAny(self.bitsize), side=side)]) + return Signature([Register('val', self.dtype, side=side)]) @staticmethod def _build_composite_state(bb: 'BloqBuilder', bits: NDArray[np.uint8]) -> Dict[str, 'SoquetT']: @@ -415,7 +419,7 @@ def _build_composite_effect( def build_composite_bloq(self, bb: 'BloqBuilder', **val: 'SoquetT') -> Dict[str, 'SoquetT']: if isinstance(self.bitsize, sympy.Expr): raise DecomposeTypeError(f'Symbolic bitsize {self.bitsize} not supported') - bits = ints_to_bits(np.array([self.val]), w=self.bitsize)[0] + bits = np.asarray(self.dtype.to_bits(self.val)) if self.state: assert not val return self._build_composite_state(bb, bits) diff --git a/qualtran/bloqs/bookkeeping/cast.py b/qualtran/bloqs/bookkeeping/cast.py index ae574953f..7b5e0a774 100644 --- a/qualtran/bloqs/bookkeeping/cast.py +++ b/qualtran/bloqs/bookkeeping/cast.py @@ -26,7 +26,6 @@ ConnectionT, DecomposeTypeError, QDType, - QFxp, Register, Side, Signature, @@ -95,12 +94,7 @@ def my_tensors( ] def on_classical_vals(self, reg: int) -> Dict[str, 'ClassicalValT']: - if isinstance(self.out_dtype, QFxp): - res = reg - elif isinstance(self.inp_dtype, QFxp): - res = int(reg) - else: - res = self.out_dtype.from_bits(self.inp_dtype.to_bits(reg)) + res = self.out_dtype.from_bits(self.inp_dtype.to_bits(reg)) return {'reg': res} def as_cirq_op(self, qubit_manager, reg: 'CirqQuregT') -> Tuple[None, Dict[str, 'CirqQuregT']]: diff --git a/qualtran/bloqs/bookkeeping/cast_test.py b/qualtran/bloqs/bookkeeping/cast_test.py index 7390258e9..ea83d9911 100644 --- a/qualtran/bloqs/bookkeeping/cast_test.py +++ b/qualtran/bloqs/bookkeeping/cast_test.py @@ -33,14 +33,14 @@ def test_cast_tensor_contraction(): def test_cast_classical_sim(): c = Cast(QInt(8), QFxp(8, 8)) (y,) = c.call_classically(reg=7) - assert y == 7 - bloq = TestCastToFrom() - (a, b) = bloq.call_classically(a=7, b=2) + assert y == 7 / 2**8 + bloq = TestCastToFrom(bitsize=8) + (a, b) = bloq.call_classically(a=7, b=QFxp(8, 8).float_to_fxp(2 / 2**8)) assert a == 7 - assert b == 9 + assert b == QFxp(8, 8).float_to_fxp(9 / 2**8) c = Cast(QFxp(8, 8), QUInt(8)) - assert c.call_classically(reg=1.2) == (1,) # type: ignore + assert c.call_classically(reg=QFxp(8, 8).float_to_fxp(1 / 2**8)) == (1,) # type: ignore def test_cast_unsiged_signed(): diff --git a/qualtran/bloqs/bookkeeping/join.py b/qualtran/bloqs/bookkeeping/join.py index e2680307c..e0c70374b 100644 --- a/qualtran/bloqs/bookkeeping/join.py +++ b/qualtran/bloqs/bookkeeping/join.py @@ -26,7 +26,6 @@ DecomposeTypeError, QBit, QDType, - QFxp, QUInt, Register, Side, @@ -34,7 +33,6 @@ ) from qualtran.bloqs.bookkeeping._bookkeeping_bloq import _BookkeepingBloq from qualtran.drawing import directional_text_box, Text, WireSymbol -from qualtran.simulation.classical_sim import bits_to_ints if TYPE_CHECKING: import quimb.tensor as qtn @@ -96,9 +94,6 @@ def my_tensors( ] def on_classical_vals(self, reg: 'NDArray[np.uint]') -> Dict[str, int]: - if isinstance(self.dtype, QFxp): - # TODO(#1095): support QFxp in classical simulation - return {'reg': bits_to_ints(reg)[0]} return {'reg': self.dtype.from_bits(reg.tolist())} def wire_symbol(self, reg: Optional[Register], idx: Tuple[int, ...] = tuple()) -> 'WireSymbol': diff --git a/qualtran/bloqs/bookkeeping/partition.py b/qualtran/bloqs/bookkeeping/partition.py index 50ce7ae6d..e04500914 100644 --- a/qualtran/bloqs/bookkeeping/partition.py +++ b/qualtran/bloqs/bookkeeping/partition.py @@ -16,6 +16,7 @@ import numpy as np from attrs import evolve, field, frozen, validators +from numpy.typing import NDArray from qualtran import ( bloq_example, @@ -24,13 +25,13 @@ ConnectionT, DecomposeTypeError, QAny, + QDType, Register, Side, Signature, ) from qualtran.bloqs.bookkeeping._bookkeeping_bloq import _BookkeepingBloq from qualtran.drawing import directional_text_box, Text, WireSymbol -from qualtran.simulation.classical_sim import bits_to_ints, ints_to_bits if TYPE_CHECKING: import quimb.tensor as qtn @@ -65,13 +66,17 @@ def __attrs_post_init__(self): if len(set(r.name for r in self.regs)) != len(self.regs): raise ValueError("Duplicate register names") + @cached_property + def lumped_dtype(self) -> QDType: + return QAny(bitsize=self.n) + @cached_property def signature(self) -> 'Signature': lumped = Side.LEFT if self.partition else Side.RIGHT partitioned = Side.RIGHT if self.partition else Side.LEFT return Signature( - [Register('x', QAny(bitsize=self.n), side=lumped)] + [Register('x', self.lumped_dtype, side=lumped)] + [evolve(reg, side=partitioned) for reg in self.regs] ) @@ -119,40 +124,35 @@ def my_tensors( def _classical_partition(self, x: 'ClassicalValT') -> Dict[str, 'ClassicalValT']: out_vals = {} - xbits = ints_to_bits(x, self.n)[0] + xbits = self.lumped_dtype.to_bits(x) start = 0 for reg in self.regs: size = int(np.prod(reg.shape + (reg.bitsize,))) bits_reg = xbits[start : start + size] if reg.shape == (): - out_vals[reg.name] = bits_to_ints(bits_reg)[0] + out_vals[reg.name] = reg.dtype.from_bits(bits_reg) else: - ints_reg = bits_to_ints( - [ - bits_reg[i * reg.bitsize : (i + 1) * reg.bitsize] - for i in range(np.prod(reg.shape)) - ] + out_vals[reg.name] = reg.dtype.from_bits_array( + np.asarray(bits_reg).reshape(reg.shape + (reg.bitsize,)) ) - out_vals[reg.name] = np.array(ints_reg).reshape(reg.shape) start += size return out_vals - def _classical_unpartition(self, **vals: 'ClassicalValT'): - out_vals = [] + def _classical_unpartition_to_bits(self, **vals: 'ClassicalValT') -> NDArray[np.uint8]: + out_vals: list[NDArray[np.uint8]] = [] for reg in self.regs: - reg_val = vals[reg.name] - if isinstance(reg_val, np.ndarray): - out_vals.append(ints_to_bits(reg_val.ravel(), reg.bitsize).ravel()) - else: - out_vals.append(ints_to_bits(reg_val, reg.bitsize)[0]) - big_int = np.concatenate(out_vals) - return {'x': bits_to_ints(big_int)[0]} + reg_val = np.asarray(vals[reg.name]) + bitstrings = reg.dtype.to_bits_array(reg_val.ravel()) + out_vals.append(bitstrings.ravel()) + return np.concatenate(out_vals) def on_classical_vals(self, **vals: 'ClassicalValT') -> Dict[str, 'ClassicalValT']: if self.partition: return self._classical_partition(vals['x']) else: - return self._classical_unpartition(**vals) + big_int_bits = self._classical_unpartition_to_bits(**vals) + big_int = self.lumped_dtype.from_bits(big_int_bits.tolist()) + return {'x': big_int} def wire_symbol(self, reg: Register, idx: Tuple[int, ...] = tuple()) -> 'WireSymbol': if reg is None: diff --git a/qualtran/bloqs/bookkeeping/split.py b/qualtran/bloqs/bookkeeping/split.py index f73a8adc8..4d124d0de 100644 --- a/qualtran/bloqs/bookkeeping/split.py +++ b/qualtran/bloqs/bookkeeping/split.py @@ -35,7 +35,6 @@ ) from qualtran.bloqs.bookkeeping._bookkeeping_bloq import _BookkeepingBloq from qualtran.drawing import directional_text_box, Text, WireSymbol -from qualtran.simulation.classical_sim import ints_to_bits if TYPE_CHECKING: import quimb.tensor as qtn @@ -88,7 +87,7 @@ def as_cirq_op(self, qubit_manager, reg: 'CirqQuregT') -> Tuple[None, Dict[str, return None, {'reg': reg.reshape((self.dtype.num_qubits, 1))} def on_classical_vals(self, reg: int) -> Dict[str, 'ClassicalValT']: - return {'reg': ints_to_bits(np.array([reg]), self.dtype.num_qubits)[0]} + return {'reg': np.asarray(self.dtype.to_bits(reg))} def my_tensors( self, incoming: Dict[str, 'ConnectionT'], outgoing: Dict[str, 'ConnectionT'] diff --git a/qualtran/bloqs/for_testing/casting.py b/qualtran/bloqs/for_testing/casting.py index 51066b1f6..d78863d36 100644 --- a/qualtran/bloqs/for_testing/casting.py +++ b/qualtran/bloqs/for_testing/casting.py @@ -34,9 +34,13 @@ @frozen class TestCastToFrom(Bloq): + bitsize: int = 4 + @cached_property def signature(self) -> Signature: - return Signature([Register('a', QUInt(4)), Register('b', QFxp(4, 4))]) + return Signature( + [Register('a', QUInt(self.bitsize)), Register('b', QFxp(self.bitsize, self.bitsize))] + ) def build_composite_bloq( self, bb: 'BloqBuilder', *, a: 'Soquet', b: 'Soquet' diff --git a/qualtran/bloqs/rotations/phase_gradient.py b/qualtran/bloqs/rotations/phase_gradient.py index 4a740115a..824c4921b 100644 --- a/qualtran/bloqs/rotations/phase_gradient.py +++ b/qualtran/bloqs/rotations/phase_gradient.py @@ -18,7 +18,6 @@ import cirq import numpy as np import sympy -from cirq._compat import cached_method from fxpmath import Fxp from numpy.typing import NDArray @@ -42,7 +41,6 @@ if TYPE_CHECKING: import quimb.tensor as qtn - from qualtran import Bloq from qualtran.resource_counting import BloqCountT, SympySymbolAllocator from qualtran.simulation.classical_sim import ClassicalValT from qualtran.symbolics import SymbolicFloat, SymbolicInt @@ -99,11 +97,15 @@ class PhaseGradientUnitary(GateWithRegisters): @cached_property def signature(self) -> 'Signature': return ( - Signature.build_from_dtypes(ctrl=QBit(), phase_grad=QFxp(self.bitsize, self.bitsize)) + Signature.build_from_dtypes(ctrl=QBit(), phase_grad=self.phase_dtype) if self.is_controlled - else Signature.build_from_dtypes(phase_grad=QFxp(self.bitsize, self.bitsize)) + else Signature.build_from_dtypes(phase_grad=self.phase_dtype) ) + @property + def phase_dtype(self) -> QFxp: + return QFxp(self.bitsize, self.bitsize) + def decompose_from_registers( self, *, context: cirq.DecompositionContext, **quregs: NDArray[cirq.Qid] # type: ignore[type-var] ) -> Iterator[cirq.OP_TREE]: @@ -123,9 +125,7 @@ def _circuit_diagram_info_(self, args: cirq.CircuitDiagramInfoArgs) -> cirq.Circ def __pow__(self, power): if power == 1: return self - return PhaseGradientUnitary( - self.bitsize, self.exponent * power, self.is_controlled, self.eps - ) + return attrs.evolve(self, exponent=self.exponent * power) def build_call_graph(self, ssa: SympySymbolAllocator) -> Set['BloqCountT']: gate = CZPowGate if self.is_controlled else ZPowGate @@ -191,9 +191,11 @@ class PhaseGradientState(GateWithRegisters): @cached_property def signature(self) -> 'Signature': - return Signature( - [Register('phase_grad', QFxp(self.bitsize, self.bitsize), side=Side.RIGHT)] - ) + return Signature([Register('phase_grad', self.phase_dtype, side=Side.RIGHT)]) + + @property + def phase_dtype(self) -> QFxp: + return QFxp(self.bitsize, self.bitsize) def decompose_from_registers( self, *, context: cirq.DecompositionContext, **quregs: NDArray[cirq.Qid] @@ -264,18 +266,15 @@ def pretty_name(self) -> str: @cached_property def signature(self) -> 'Signature': return ( - Signature.build_from_dtypes( - ctrl=QBit(), - x=QFxp(self.x_bitsize, self.x_bitsize, signed=False), - phase_grad=QFxp(self.phase_bitsize, self.phase_bitsize, signed=False), - ) + Signature.build_from_dtypes(ctrl=QBit(), x=self.x_dtype, phase_grad=self.phase_dtype) if self.controlled_by is not None - else Signature.build_from_dtypes( - x=QFxp(self.x_bitsize, self.x_bitsize, signed=False), - phase_grad=QFxp(self.phase_bitsize, self.phase_bitsize, signed=False), - ) + else Signature.build_from_dtypes(x=self.x_dtype, phase_grad=self.phase_dtype) ) + @cached_property + def x_dtype(self) -> QFxp: + return QFxp(self.x_bitsize, self.x_bitsize, signed=False) + @cached_property def phase_dtype(self) -> QFxp: return QFxp(self.phase_bitsize, self.phase_bitsize, signed=False) @@ -292,37 +291,52 @@ def registers(self) -> Sequence[Union[int, Sequence[int]]]: def with_registers(self, *new_registers: Union[int, Sequence[int]]): raise NotImplementedError("not needed.") - @cached_method - def scaled_val(self, x: int) -> int: - """Computes `phase_grad + x` using fixed point arithmetic.""" - x_width = self.x_bitsize + self.right_shift - x_fxp = _fxp(x / 2**x_width, x_width).like(_fxp(0, self.phase_bitsize)).astype(float) - return int(x_fxp.astype(float) * 2**self.phase_bitsize) - def apply(self, *args) -> Tuple[Union[int, np.integer, NDArray[np.integer]], ...]: if self.controlled_by is not None: ctrl, x, phase_grad = args - out = self.on_classical_vals(ctrl=ctrl, x=x, phase_grad=phase_grad) - return out['ctrl'], out['x'], out['phase_grad'] + else: + ctrl = None + x, phase_grad = args + + if self.controlled_by == ctrl: + x_fxp = self.x_dtype.float_to_fxp(x, raw=True) + phase_grad_fxp = self.phase_dtype.float_to_fxp(phase_grad, raw=True) - x, phase_grad = args - out = self.on_classical_vals(x=x, phase_grad=phase_grad) - return out['x'], out['phase_grad'] + ctrl_kwarg = {} if self.controlled_by is None else {'ctrl': ctrl} + out = self.on_classical_vals(x=x_fxp, phase_grad=phase_grad_fxp, **ctrl_kwarg) + + phase_grad_out_fxp = out['phase_grad'] + assert isinstance(phase_grad_out_fxp, Fxp) + phase_grad_out = int(phase_grad_out_fxp.raw()) + else: + phase_grad_out = phase_grad + + if ctrl is not None: + return ctrl, x, phase_grad_out + else: + return x, phase_grad_out def on_classical_vals(self, **kwargs) -> Dict[str, 'ClassicalValT']: - x, phase_grad = kwargs['x'], kwargs['phase_grad'] - if self.controlled_by is not None: - ctrl = kwargs['ctrl'] - if ctrl == self.controlled_by: - phase_grad_out = (phase_grad + self.sign * self.scaled_val(x)) % ( - 2**self.phase_bitsize - ) + x: Fxp = kwargs['x'] + phase_grad: Fxp = kwargs['phase_grad'] + + if self.controlled_by is None or self.controlled_by == kwargs['ctrl']: + # widen appropriately so that right shifting does not drop necessary bits + x = x.like( + QFxp( + x.n_int + phase_grad.n_frac, phase_grad.n_frac, phase_grad.signed + ).fxp_dtype_template + ) + scaled_x = x >> self.right_shift + if self.sign == 1: + phase_grad_out = phase_grad + scaled_x else: - phase_grad_out = phase_grad - return {'ctrl': ctrl, 'x': x, 'phase_grad': phase_grad_out} + phase_grad_out = phase_grad - scaled_x + else: + phase_grad_out = phase_grad.copy() - phase_grad_out = (phase_grad + self.sign * self.scaled_val(x)) % (2**self.phase_bitsize) - return {'x': x, 'phase_grad': phase_grad_out} + kwargs['phase_grad'] = phase_grad_out + return kwargs def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']: num_toffoli = self.phase_bitsize - 2 @@ -331,21 +345,8 @@ def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']: return {(Toffoli(), num_toffoli)} - def adjoint(self) -> 'Bloq': - return AddIntoPhaseGrad( - self.x_bitsize, - self.phase_bitsize, - self.right_shift, - sign=-1 * self.sign, - controlled_by=self.controlled_by, - ) - - def __pow__(self, power): - if power == 1: - return self - if power == -1: - return self.adjoint() - raise NotImplementedError("AddIntoPhaseGrad.__pow__ defined only for powers +1/-1.") + def adjoint(self) -> 'AddIntoPhaseGrad': + return attrs.evolve(self, sign=-self.sign) def my_tensors( self, incoming: Dict[str, 'ConnectionT'], outgoing: Dict[str, 'ConnectionT'] @@ -365,37 +366,30 @@ def _add_into_phase_grad() -> AddIntoPhaseGrad: def _fxp(x: float, n: 'SymbolicInt') -> Fxp: - """When 0 <= x < 1, constructs an n-bit fixed point representation with nice properties. - - Specifically, - - op_sizing='same' and const_op_sizing='same' ensure that the returned object is not resized - to a bigger fixed point number when doing operations with other Fxp objects. - - shifting='trunc' ensures that when shifting the Fxp integer to left / right; the digits are - truncated and no rounding occurs - - overflow='wrap' ensures that when performing operations where result overflows, the overflowed - digits are simply discarded. - """ + """When 0 <= x < 1, constructs an n-bit fixed point representation with nice properties.""" assert 0 <= x < 1 - return Fxp( - x, - dtype=f'fxp-u{n}/{n}', - op_sizing='same', - const_op_sizing='same', - shifting='trunc', - overflow='wrap', - ) + return QFxp(n, n).float_to_fxp(x, require_exact=False) def _mul_via_repeated_add(x_fxp: Fxp, gamma_fxp: Fxp, out: int) -> Fxp: """Multiplication via repeated additions algorithm described in Appendix D5""" res = _fxp(0, out) + x_fxp = x_fxp.copy() + x_fxp.resize(n_int=x_fxp.n_int, n_frac=max(out, x_fxp.n_frac)) for i, bit in enumerate(gamma_fxp.bin()): if bit == '0': continue shift = gamma_fxp.n_int - i - 1 # Left/Right shift by `shift` bits. - res += x_fxp << shift if shift > 0 else x_fxp >> abs(shift) + # Issue: Fxp << and >> does not preserve config, so we use `* 2**shift` instead. + x_shifted = x_fxp.copy() + if shift > 0: + x_shifted.set_val(x_fxp.val << np.array(shift, dtype=x_fxp.val.dtype), raw=True) + else: + x_shifted.set_val(x_fxp.val >> np.array(-shift, dtype=x_fxp.val.dtype), raw=True) + assert x_shifted.overflow == 'wrap' and x_shifted.shifting == 'trunc' + res += x_shifted return res @@ -459,34 +453,10 @@ def phase_dtype(self) -> QFxp: return QFxp(self.phase_bitsize, self.phase_bitsize, signed=False) @cached_property - def gamma_fxp(self) -> Fxp: - return Fxp(abs(self.gamma), dtype=self.gamma_dtype.fxp_dtype_str) - - @cached_method - def scaled_val(self, x: int) -> int: - """Computes `x*self.gamma` using fixed point arithmetic.""" - if isinstance(self.gamma, sympy.Expr): - raise ValueError(f'Symbolic gamma {self.gamma} not allowed') - if isinstance(self.phase_dtype.bitsize, sympy.Expr): - raise ValueError(f'Symbolic phase bitsize {self.phase_dtype.bitsize} not allowed') - sign = np.sign(self.gamma) - # `x` is an integer because we currently represent each bitstring as an integer during simulations. - # However, `x` should be interpreted as per the fixed point specification given in self.x_dtype. - # If `self.x_dtype` uses `n_frac` bits to represent the fractional part, `x` should be divided by - # 2**n_frac (in other words, right shifted by n_frac) - x_fxp = Fxp(x / 2**self.x_dtype.num_frac, dtype=self.x_dtype.fxp_dtype_str) - # Similarly, `self.gamma` should be represented as a fixed point number using appropriate number - # of bits for integer and fractional part. This is done in self.gamma_fxp - # Compute the result = x_fxp * gamma_fxp - result = _mul_via_repeated_add(x_fxp, self.gamma_fxp, self.phase_dtype.bitsize) - # Since the phase gradient register is a fixed point register with `n_word=0`, we discard the integer - # part of `result`. This is okay because the adding `x.y` to the phase gradient register will impart - # a phase of `exp(i * 2 * np.pi * x.y)` which is same as `exp(i * 2 * np.pi * y)` - assert 0 <= result < 1 and result.dtype == self.phase_dtype.fxp_dtype_str - # Convert the `self.phase_bitsize`-bit fraction into back to an integer and return the result. - # Sign of `gamma` affects whether we add or subtract into the phase gradient register and thus - # can be ignored during the fixed point arithmetic analysis. - return int(np.floor(result.astype(float) * 2**self.phase_dtype.bitsize) * sign) + def abs_gamma_fxp(self) -> Fxp: + if is_symbolic(self.gamma): + raise ValueError(f"Cannot compute Fxp for symbolic {self.gamma=}") + return self.gamma_dtype.float_to_fxp(abs(self.gamma), require_exact=False) def decompose_from_registers( self, *, context: cirq.DecompositionContext, **quregs: NDArray[cirq.Qid] @@ -495,16 +465,16 @@ def decompose_from_registers( raise ValueError(f'Symbolic gamma {self.gamma} not allowed') x, phase_grad = quregs['x'], quregs['phase_grad'] sign = int(np.sign(self.gamma)) - for i, bit in enumerate(self.gamma_fxp.bin()): + for i, bit in enumerate(self.abs_gamma_fxp.bin()): if bit == '0': continue - shift = self.gamma_fxp.n_int - i - 1 + shift = self.abs_gamma_fxp.n_int - i - 1 if 0 <= shift < self.x_dtype.num_frac: # Left shift by `shift` bits / multiply by 2**shift yield AddIntoPhaseGrad( self.x_dtype.num_frac - shift, self.phase_bitsize, sign=sign ).on_registers(x=x[shift + self.x_dtype.num_int :], phase_grad=phase_grad) - elif -len(phase_grad) - self.x_dtype.num_int < shift < 0: + elif -self.phase_dtype.bitsize - self.x_dtype.num_int < shift < 0: # Right shift by `shift` bits / divide by 2**shift shift = abs(shift) x_bitsize = self.x_dtype.num_frac + min(shift, self.x_dtype.num_int) @@ -520,21 +490,37 @@ def apply( ) -> Tuple[ Union[int, np.integer, NDArray[np.integer]], Union[int, np.integer, NDArray[np.integer]] ]: - out = self.on_classical_vals(x=x, phase_grad=phase_grad) - return out['x'], out['phase_grad'] + x_fxp = self.x_dtype.float_to_fxp(x, raw=True) + phase_grad_fxp = self.phase_dtype.float_to_fxp(phase_grad, raw=True) - def on_classical_vals(self, x: int, phase_grad: int) -> Dict[str, 'ClassicalValT']: - phase_grad_out = (phase_grad + self.scaled_val(x)) % 2**self.phase_bitsize + out = self.on_classical_vals(x=x_fxp, phase_grad=phase_grad_fxp) + phase_grad_out = int(out['phase_grad'].raw()) + + return x, phase_grad_out + + def on_classical_vals(self, x: Fxp, phase_grad: Fxp) -> Dict[str, Fxp]: + if is_symbolic(self.phase_bitsize): + raise ValueError(f"Cannot classically simulate with symbolic {self.phase_bitsize=}") + if is_symbolic(self.gamma): + raise ValueError(f"Cannot classically simulate with symbolic {self.gamma=}") + + scaled_x = _mul_via_repeated_add( + x.like(phase_grad), gamma_fxp=self.abs_gamma_fxp, out=self.phase_bitsize + ) + if np.sign(self.gamma) == -1: + phase_grad_out = phase_grad - scaled_x + else: + phase_grad_out = phase_grad + scaled_x return {'x': x, 'phase_grad': phase_grad_out} def build_call_graph(self, ssa: 'SympySymbolAllocator') -> Set['BloqCountT']: num_additions = (self.gamma_dtype.bitsize + 2) // 2 if not isinstance(self.gamma, sympy.Basic): num_additions_naive = 0 - for i, bit in enumerate(self.gamma_fxp.bin()): + for i, bit in enumerate(self.abs_gamma_fxp.bin()): if bit == '0': continue - shift = self.gamma_fxp.n_int - i - 1 + shift = self.abs_gamma_fxp.n_int - i - 1 if -(self.phase_bitsize + self.x_dtype.num_int) < shift < self.x_dtype.num_frac: num_additions_naive += 1 num_additions = min(num_additions_naive, num_additions) diff --git a/qualtran/bloqs/rotations/phase_gradient_test.py b/qualtran/bloqs/rotations/phase_gradient_test.py index bb02845ee..c1f9386d5 100644 --- a/qualtran/bloqs/rotations/phase_gradient_test.py +++ b/qualtran/bloqs/rotations/phase_gradient_test.py @@ -77,38 +77,88 @@ def test_phase_gradient_gate(n: int, exponent, controlled): assert np.allclose(cirq.unitary(bloq), cirq.unitary(cirq_gate), atol=eps) -def test_add_into_phase_grad(): +@pytest.mark.slow +def test_add_into_phase_grad_classical_sim(): + from qualtran.bloqs.rotations.phase_gradient import _fxp + + x_bit, phase_bit = 4, 7 + bloq = AddIntoPhaseGrad(x_bit, phase_bit) + for x in range(2**x_bit): + for phase_grad in range(2**phase_bit): + phase_fxp = _fxp(phase_grad / 2**phase_bit, phase_bit) + x_fxp = _fxp(x / 2**x_bit, x_bit) + assert bloq.call_classically(x=x_fxp, phase_grad=phase_fxp) == ( + x_fxp, + phase_fxp + x_fxp, + ) + + +@pytest.mark.slow +def test_add_into_phase_grad_unitary(): from qualtran.bloqs.rotations.phase_gradient import _fxp x_bit, phase_bit = 4, 7 bloq = AddIntoPhaseGrad(x_bit, phase_bit) + + # Prepare basis states mapping for cirq-style simulation. basis_map: Dict[int, int] = {} for x in range(2**x_bit): for phase_grad in range(2**phase_bit): phase_fxp = _fxp(phase_grad / 2**phase_bit, phase_bit) - x_fxp = _fxp(x / 2**x_bit, x_bit).like(phase_fxp) + x_fxp = _fxp(x / 2**x_bit, x_bit) phase_grad_out = int((phase_fxp + x_fxp).astype(float) * 2**phase_bit) - # Test Bloq style classical simulation. - assert bloq.call_classically(x=x, phase_grad=phase_grad) == (x, phase_grad_out) - # Prepare basis states mapping for cirq-style simulation. + input_state = int(f'{x:0{x_bit}b}' + f'{phase_grad:0{phase_bit}b}', 2) output_state = int(f'{x:0{x_bit}b}' + f'{phase_grad_out:0{phase_bit}b}', 2) basis_map[input_state] = output_state + # Test cirq style simulation. num_bits = x_bit + phase_bit assert len(basis_map) == len(set(basis_map.values())) circuit = cirq.Circuit(bloq.on(*cirq.LineQubit.range(num_bits))) cirq.testing.assert_equivalent_computational_basis_map(basis_map, circuit) + + +def test_add_into_phase_grad_t_complexity(): + x_bit, phase_bit = 4, 7 + bloq = AddIntoPhaseGrad(x_bit, phase_bit) ((toffoli, n),) = bloq.bloq_counts().items() assert bloq.t_complexity() == n * toffoli.t_complexity() +@pytest.mark.slow +@pytest.mark.parametrize('controlled', [0, 1]) +def test_add_into_phase_grad_controlled_classical_sim(controlled: int): + from qualtran.bloqs.rotations.phase_gradient import _fxp + + x_bit, phase_bit = 4, 7 + bloq = AddIntoPhaseGrad(x_bit, phase_bit, controlled_by=controlled) + + for control in range(2): + for x in range(2**x_bit): + for phase_grad in range(2**phase_bit): + phase_fxp = _fxp(phase_grad / 2**phase_bit, phase_bit) + x_fxp = _fxp(x / 2**x_bit, x_bit) + if control == controlled: + phase_grad_out_fxp = phase_fxp + x_fxp + else: + phase_grad_out_fxp = phase_fxp + assert bloq.call_classically(ctrl=control, x=x_fxp, phase_grad=phase_fxp) == ( + control, + x_fxp, + phase_grad_out_fxp, + ) + + +@pytest.mark.slow @pytest.mark.parametrize('controlled', [0, 1]) def test_add_into_phase_grad_controlled(controlled: int): from qualtran.bloqs.rotations.phase_gradient import _fxp x_bit, phase_bit = 4, 7 bloq = AddIntoPhaseGrad(x_bit, phase_bit, controlled_by=controlled) + + # Prepare basis states mapping for cirq-style simulation. basis_map: Dict[int, int] = {} num_bits = 1 + x_bit + phase_bit expected_unitary = np.zeros((2**num_bits, 2**num_bits)) @@ -116,18 +166,12 @@ def test_add_into_phase_grad_controlled(controlled: int): for x in range(2**x_bit): for phase_grad in range(2**phase_bit): phase_fxp = _fxp(phase_grad / 2**phase_bit, phase_bit) - x_fxp = _fxp(x / 2**x_bit, x_bit).like(phase_fxp) + x_fxp = _fxp(x / 2**x_bit, x_bit) if control == controlled: - phase_grad_out = int((phase_fxp + x_fxp).astype(float) * 2**phase_bit) + phase_grad_out_fxp = phase_fxp + x_fxp + phase_grad_out = int(phase_grad_out_fxp.astype(float) * 2**phase_bit) else: phase_grad_out = phase_grad - # Test Bloq style classical simulation. - assert bloq.call_classically(ctrl=control, x=x, phase_grad=phase_grad) == ( - control, - x, - phase_grad_out, - ) - # Prepare basis states mapping for cirq-style simulation. input_state = int( f'{control}' + f'{x:0{x_bit}b}' + f'{phase_grad:0{phase_bit}b}', 2 ) @@ -136,36 +180,51 @@ def test_add_into_phase_grad_controlled(controlled: int): ) basis_map[input_state] = output_state expected_unitary[output_state, input_state] = 1 + # Test cirq style simulation. assert len(basis_map) == len(set(basis_map.values())) circuit = cirq.Circuit(bloq.on(*cirq.LineQubit.range(num_bits))) np.testing.assert_allclose(circuit.unitary(), expected_unitary, atol=1e-8) +_ADD_SCALED_VAL_INTO_PHASE_REG_EXAMPLES: list[AddScaledValIntoPhaseReg] = [ + AddScaledValIntoPhaseReg.from_bitsize(4, 7, 0.123, 6), + AddScaledValIntoPhaseReg.from_bitsize(2, 8, 1.3868682, 8), + AddScaledValIntoPhaseReg.from_bitsize(4, 9, -9.0949456, 5), + AddScaledValIntoPhaseReg.from_bitsize(6, 4, 2.5, 2), + AddScaledValIntoPhaseReg(QFxp(4, 4, signed=False), 4, 1.3868682, QFxp(8, 7, signed=False)), +] + + @pytest.mark.slow -@pytest.mark.parametrize( - 'bloq', - [ - AddScaledValIntoPhaseReg.from_bitsize(4, 7, 0.123, 6), - AddScaledValIntoPhaseReg.from_bitsize(2, 8, 1.3868682, 8), - AddScaledValIntoPhaseReg.from_bitsize(4, 9, -19.0949456, 5), - AddScaledValIntoPhaseReg.from_bitsize(6, 4, 2.5, 2), - AddScaledValIntoPhaseReg(QFxp(4, 0, signed=False), 4, 1.3868682, QFxp(8, 7, signed=False)), - ], -) -def test_add_scaled_val_into_phase_reg(bloq): +@pytest.mark.parametrize('bloq', _ADD_SCALED_VAL_INTO_PHASE_REG_EXAMPLES) +def test_add_scaled_val_into_phase_reg_classical_sim(bloq): cbloq = bloq.decompose_bloq() - for x in range(2**bloq.x_dtype.bitsize): - for phase_grad in range(2**bloq.phase_bitsize): + for x in bloq.x_dtype.get_classical_domain(): + for phase_grad in bloq.phase_dtype.get_classical_domain(): d = {'x': x, 'phase_grad': phase_grad} c1 = bloq.on_classical_vals(**d) c2 = cbloq.on_classical_vals(**d) assert c1 == c2, f'{d=}, {c1=}, {c2=}' + + +@pytest.mark.parametrize( + 'bloq', + [ + pytest.param(bloq, marks=pytest.mark.slow if bloq.num_qubits() > 11 else ()) + for bloq in _ADD_SCALED_VAL_INTO_PHASE_REG_EXAMPLES + ], +) +def test_add_scaled_val_into_phase_reg_unitary(bloq: AddScaledValIntoPhaseReg): bloq_unitary = cirq.unitary(bloq) op = GateHelper(bloq).operation circuit = cirq.Circuit(cirq.I.on_each(*op.qubits), cirq.decompose_once(op)) decomposed_unitary = circuit.unitary(qubit_order=op.qubits) np.testing.assert_allclose(bloq_unitary, decomposed_unitary) + + +@pytest.mark.parametrize('bloq', _ADD_SCALED_VAL_INTO_PHASE_REG_EXAMPLES) +def test_add_scaled_val_into_phase_reg_unitary_t_complexity(bloq: AddScaledValIntoPhaseReg): ((add_into_phase, n),) = bloq.bloq_counts().items() assert bloq.t_complexity() == n * add_into_phase.t_complexity() diff --git a/qualtran/bloqs/rotations/phasing_via_cost_function_test.py b/qualtran/bloqs/rotations/phasing_via_cost_function_test.py index ebb8bdc18..e7fa20dad 100644 --- a/qualtran/bloqs/rotations/phasing_via_cost_function_test.py +++ b/qualtran/bloqs/rotations/phasing_via_cost_function_test.py @@ -127,7 +127,6 @@ def test_hamming_weight_phasing_using_phase_via_cost_function( @attrs.frozen class TestSquarePhasing(GateWithRegisters): bitsize: int - normalize: bool gamma: float eps: float use_phase_gradient: bool = False @@ -138,9 +137,7 @@ def signature(self) -> 'Signature': @property def cost_reg(self) -> Register: - return Register( - 'result', QFxp(2 * self.bitsize, 2 * self.bitsize * int(self.normalize), signed=False) - ) + return Register('result', QFxp(2 * self.bitsize, 2 * self.bitsize, signed=False)) @property def phase_gradient_oracle(self) -> QvrPhaseGradient: @@ -169,30 +166,17 @@ def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: 'SoquetT') -> Dict[str return soqs -@pytest.mark.slow -@pytest.mark.parametrize('normalize_cost_function', [True, False]) -@pytest.mark.parametrize('use_phase_gradient', [True, False]) -@pytest.mark.parametrize('gamma, eps', [(0.1, 5e-2), (1.20345, 5e-2), (-1.1934341, 5e-2)]) +@pytest.mark.parametrize('use_phase_gradient', [pytest.param(True, marks=pytest.mark.slow), False]) +@pytest.mark.parametrize('gamma, eps', [(0.125, 5e-2), (1.1875, 5e-2), (-1.375, 5e-2)]) @pytest.mark.parametrize('n', [2]) def test_square_phasing_via_phase_gradient( - n: int, gamma: float, eps: float, use_phase_gradient: bool, normalize_cost_function: bool + n: int, gamma: float, eps: float, use_phase_gradient: bool ): initial_state = np.array([1 / np.sqrt(2**n)] * 2**n) - normalization_factor = 1 if normalize_cost_function else 4**n - phases = np.array( - [ - np.exp(1j * 2 * np.pi * gamma * x**2 * normalization_factor / 4**n) - for x in range(2**n) - ] - ) + phases = np.array([np.exp(1j * 2 * np.pi * gamma * x**2 / 4**n) for x in range(2**n)]) expected_final_state = np.multiply(initial_state, phases) - test_bloq_one = TestSquarePhasing( - n, True, gamma * normalization_factor, eps, use_phase_gradient - ) - test_bloq_two = TestSquarePhasing( - n, False, gamma * normalization_factor / (4**n), eps, use_phase_gradient - ) - for test_bloq in [test_bloq_one, test_bloq_two]: + test_bloq = TestSquarePhasing(n, gamma, eps, use_phase_gradient) + for test_bloq in [test_bloq]: bb = BloqBuilder() a = bb.allocate(n) a = bb.add(OnEach(n, Hadamard()), q=a) diff --git a/qualtran/bloqs/rotations/quantum_variable_rotation.py b/qualtran/bloqs/rotations/quantum_variable_rotation.py index 1559b7e94..73fb7ebfb 100644 --- a/qualtran/bloqs/rotations/quantum_variable_rotation.py +++ b/qualtran/bloqs/rotations/quantum_variable_rotation.py @@ -235,12 +235,13 @@ def find_optimal_phase_grad_size(gamma_fxp: Fxp, cost_dtype: QFxp, eps: float) - from qualtran.bloqs.rotations.phase_gradient import _mul_via_repeated_add cost_val = (2**cost_dtype.bitsize - 1) / (2**cost_dtype.num_frac) - cost_fxp = Fxp(cost_val, dtype=cost_dtype.fxp_dtype_str) + cost_fxp = cost_dtype.float_to_fxp(cost_val) expected_val = (gamma_fxp.get_val() * cost_val) % 1 def is_good_phase_grad_size(phase_bitsize: int): res = _mul_via_repeated_add(cost_fxp, gamma_fxp, phase_bitsize) - return np.abs(res.get_val() - expected_val) <= eps + actual_val = res.get_val() % 1 + return np.abs(actual_val - expected_val) <= eps low, high, ans = 1, 100, None while low <= high: @@ -461,7 +462,10 @@ def num_additions(self) -> int: @cached_property def gamma_fxp(self) -> Fxp: - return Fxp(abs(self.gamma), dtype=self.gamma_dtype.fxp_dtype_str) + if is_symbolic(self.gamma): + raise ValueError(f"Cannot compute Fxp from symbolic {self.gamma=}") + + return self.gamma_dtype.float_to_fxp(abs(self.gamma), require_exact=False) @cached_property def gamma_dtype(self) -> QFxp: @@ -472,7 +476,7 @@ def gamma_dtype(self) -> QFxp: # The reference assumes that cost register always stores a fraction between [0, 1). We # do not have this assumption and therefore, we also need to add self.cost_dtype.num_int # to the gamma bitsize. - n_int = smax(0, bit_length(sympy.Abs(self.gamma))) + n_int = smax(0, bit_length(sabs(self.gamma))) n_frac = self.cost_dtype.num_int + self.b_phase return QFxp(bitsize=n_int + n_frac, num_frac=n_frac, signed=False) diff --git a/qualtran/simulation/classical_sim.py b/qualtran/simulation/classical_sim.py index 63af89cca..bd00b6fa0 100644 --- a/qualtran/simulation/classical_sim.py +++ b/qualtran/simulation/classical_sim.py @@ -19,6 +19,7 @@ import networkx as nx import numpy as np import sympy +from fxpmath import Fxp from numpy.typing import NDArray from qualtran import ( @@ -34,35 +35,7 @@ ) from qualtran._infra.composite_bloq import _binst_to_cxns -ClassicalValT = Union[int, np.integer, NDArray[np.integer]] - - -def bits_to_ints(bitstrings: Union[Sequence[int], NDArray[np.uint]]) -> NDArray[np.uint]: - """Returns the integer specified by the given big-endian bitstrings. - - Args: - bitstrings: A bitstring or array of bitstrings, each of which has the 1s bit (LSB) at the end. - Returns: - An array of integers; one for each bitstring. - """ - bitstrings = np.atleast_2d(bitstrings) - if bitstrings.shape[1] > 64: - raise NotImplementedError() - basis = 2 ** np.arange(bitstrings.shape[1] - 1, 0 - 1, -1, dtype=np.uint64) - return np.sum(basis * bitstrings, axis=1, dtype=np.uint64) - - -def ints_to_bits( - x: Union[int, np.integer, Sequence[int], NDArray[np.integer]], w: int -) -> NDArray[np.uint8]: - """Returns the big-endian bitstrings specified by the given integers. - - Args: - x: An integer or array of unsigned integers. - w: The bit width of the returned bitstrings. - """ - x = np.atleast_1d(x) - return np.array([list(map(int, np.binary_repr(v, width=w))) for v in x], dtype=np.uint8) +ClassicalValT = Union[int, np.integer, NDArray[np.integer], Fxp] def _get_in_vals( diff --git a/qualtran/simulation/classical_sim_test.py b/qualtran/simulation/classical_sim_test.py index 398d640e3..5091dbe12 100644 --- a/qualtran/simulation/classical_sim_test.py +++ b/qualtran/simulation/classical_sim_test.py @@ -15,7 +15,6 @@ import itertools from typing import Dict -import cirq import numpy as np import pytest from attrs import frozen @@ -26,53 +25,12 @@ from qualtran.simulation.classical_sim import ( _update_assign_from_vals, add_ints, - bits_to_ints, call_cbloq_classically, ClassicalValT, - ints_to_bits, ) from qualtran.testing import execute_notebook -def test_bits_to_int(): - rs = np.random.RandomState(52) - bitstrings = rs.choice([0, 1], size=(100, 23)) - - nums = bits_to_ints(bitstrings) - assert nums.dtype == np.uint64 - assert nums.shape == (100,) - - for num, bs in zip(nums, bitstrings): - ref_num = cirq.big_endian_bits_to_int(bs.tolist()) - assert num == ref_num - - # check one input bitstring instead of array of input bitstrings. - (num,) = bits_to_ints([1, 0]) - assert num == 2 - - -def test_int_to_bits(): - rs = np.random.RandomState(52) - nums = rs.randint(0, 2**23 - 1, size=(100,), dtype=np.uint64) - bitstrings = ints_to_bits(nums, w=23) - assert bitstrings.shape == (100, 23) - - nums = rs.randint(-(2**22), 2**22, size=(100,), dtype=np.int64) - bitstrings = ints_to_bits(nums, w=23) - assert bitstrings.shape == (100, 23) - - for num, bs in zip(nums, bitstrings): - ref_bs = cirq.big_endian_int_to_bits(int(num), bit_count=23) - np.testing.assert_array_equal(ref_bs, bs) - - # check one input int - (bitstring,) = ints_to_bits(2, w=8) - assert bitstring.tolist() == [0, 0, 0, 0, 0, 0, 1, 0] - - bitstring = ints_to_bits([31, -1], w=6) - assert bitstring.tolist() == [[0, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1]] - - def test_dtype_validation(): # set up mocks for `_update_assign_from_vals` soq_assign: Dict[Soquet, ClassicalValT] = {} # gets assigned to; we discard in this test.