diff --git a/recirq/fermi_hubbard/parameters.py b/recirq/fermi_hubbard/parameters.py
index 29f4131b..f0618c4c 100644
--- a/recirq/fermi_hubbard/parameters.py
+++ b/recirq/fermi_hubbard/parameters.py
@@ -18,7 +18,7 @@
 
 import abc
 from itertools import product
-from numbers import Number
+from numbers import Number, Real
 
 import cirq
 import numpy as np
@@ -30,8 +30,6 @@
     ZigZagLayout
 )
 
-Real = Union[int, float]
-
 
 @dataclass(init=False)
 class Hamiltonian:
@@ -102,45 +100,45 @@ def interactions_count(self) -> int:
 
     @property
     def j_array(self) -> np.ndarray:
-        if isinstance(self.j, tuple):
-            return np.array(self.j)
-        return np.full(self.interactions_count, self.j)
+        if isinstance(self.j, Real):
+            return np.full(self.interactions_count, self.j)
+        return np.array(self.j)
 
     @property
     def u_array(self) -> np.ndarray:
-        if isinstance(self.u, tuple):
-            return np.array(self.u)
-        return np.full(self.sites_count, self.u)
+        if isinstance(self.u, Real):
+            return np.full(self.sites_count, self.u)
+        return np.array(self.u)
 
     @property
     def v_array(self) -> np.ndarray:
-        if isinstance(self.v, tuple):
-            return np.array(self.v)
-        return np.full(self.interactions_count, self.v)
+        if isinstance(self.v, Real):
+            return np.full(self.interactions_count, self.v)
+        return np.array(self.v)
 
     @property
     def local_charge_array(self) -> np.ndarray:
-        if isinstance(self.local_charge, tuple):
-            return np.array(self.local_charge)
-        return np.full(self.sites_count, self.local_charge)
+        if isinstance(self.local_charge, Real):
+            return np.full(self.sites_count, self.local_charge)
+        return np.array(self.local_charge)
 
     @property
     def local_spin_array(self) -> np.ndarray:
-        if isinstance(self.local_spin, tuple):
-            return np.array(self.local_spin)
-        return np.full(self.sites_count, self.local_spin)
+        if isinstance(self.local_spin, Real):
+            return np.full(self.sites_count, self.local_spin)
+        return np.array(self.local_spin)
 
     @property
     def mu_up_array(self) -> np.ndarray:
-        if isinstance(self.mu_up, tuple):
-            return np.array(self.mu_up)
-        return np.full(self.sites_count, self.mu_up)
+        if isinstance(self.mu_up, Real):
+            return np.full(self.sites_count, self.mu_up)
+        return np.array(self.mu_up)
 
     @property
     def mu_down_array(self) -> np.ndarray:
-        if isinstance(self.mu_down, tuple):
-            return np.array(self.mu_down)
-        return np.full(self.sites_count, self.mu_down)
+        if isinstance(self.mu_down, Real):
+            return np.full(self.sites_count, self.mu_down)
+        return np.array(self.mu_down)
 
     @property
     def local_up_array(self):
@@ -257,7 +255,7 @@ def get_amplitudes(self, sites_count: int) -> np.ndarray:
         if sites_count != len(self.amplitudes):
             raise ValueError(f'Fixed single particle not compatible with '
                              f'{sites_count} sites')
-        return np.array(self.potential)
+        return np.array(self.amplitudes)
 
     def _json_dict_(self):
         return cirq.dataclass_json_dict(self)
@@ -638,10 +636,10 @@ def _potential_to_quadratic_hamiltonian(
 ) -> openfermion.QuadraticHamiltonian:
     sites_count = len(potential)
 
-    if isinstance(j, Iterable):
-        j = np.array(j)
-    else:
+    if isinstance(j, Real):
         j = np.full(sites_count - 1, j)
+    else:
+        j = np.array(j)
 
     if len(j) != sites_count - 1:
         raise ValueError('Hopping coefficient size incompatible with potential')
diff --git a/recirq/kpz/experiment.py b/recirq/kpz/experiment.py
index 50f197f2..8edb8d88 100644
--- a/recirq/kpz/experiment.py
+++ b/recirq/kpz/experiment.py
@@ -45,6 +45,7 @@
 """
 
 from typing import Iterator, List, Union, Optional
+from numbers import Real
 
 import cirq
 import numpy as np
@@ -57,7 +58,7 @@
 rng = np.random.default_rng()
 
 
-def _dec_to_binary_right(d: Union[np.ndarray, int], n: int) -> Union[np.ndarray, int]:
+def _dec_to_binary_right(d: Union[np.ndarray, Real], n: int) -> Union[np.ndarray, Real]:
     i = np.arange(n // 2)
     return (
         np.floor(np.outer(d, 1 / 2**i)) - np.floor(np.outer(d, 1 / 2 ** (i + 1))) * 2
@@ -103,25 +104,25 @@ def __init__(self, prob_right: np.ndarray, initial_states: np.ndarray):
         self.skewness = self._skewness()
         self.kurtosis = self._kurtosis()
 
-    def _mean(self) -> float:
+    def _mean(self) -> Real:
         return (
             self.transferred_magnetization_probs @ self.transferred_magnetization_vals
         )
 
-    def _variance(self) -> float:
+    def _variance(self) -> Real:
         return (
             self.transferred_magnetization_probs
             @ (self.transferred_magnetization_vals - self.mean) ** 2
         )
 
-    def _skewness(self) -> float:
+    def _skewness(self) -> Real:
         return (
             self.transferred_magnetization_probs
             @ (self.transferred_magnetization_vals - self.mean) ** 3
             / self.variance ** (3 / 2)
         )
 
-    def _kurtosis(self) -> float:
+    def _kurtosis(self) -> Real:
         return (
             self.transferred_magnetization_probs
             @ (self.transferred_magnetization_vals - self.mean) ** 4
@@ -129,20 +130,20 @@ def _kurtosis(self) -> float:
             - 3
         )
 
-    def _mean_excluding_i(self, i: int) -> float:
+    def _mean_excluding_i(self, i: int) -> Real:
         p = np.mean(
             np.delete(self.transferred_magnetization_probs_all, i, axis=0), axis=0
         )
         return p @ self.transferred_magnetization_vals
 
-    def _variance_excluding_i(self, i: int) -> float:
+    def _variance_excluding_i(self, i: int) -> Real:
         p = np.mean(
             np.delete(self.transferred_magnetization_probs_all, i, axis=0), axis=0
         )
         mean_i = p @ self.transferred_magnetization_vals
         return p @ (self.transferred_magnetization_vals - mean_i) ** 2
 
-    def _skew_excluding_i(self, i: int) -> float:
+    def _skew_excluding_i(self, i: int) -> Real:
         p = np.mean(
             np.delete(self.transferred_magnetization_probs_all, i, axis=0), axis=0
         )
@@ -154,7 +155,7 @@ def _skew_excluding_i(self, i: int) -> float:
             / variance_i ** (3 / 2)
         )
 
-    def _kurtosis_excluding_i(self, i: int) -> float:
+    def _kurtosis_excluding_i(self, i: int) -> Real:
         p = np.mean(
             np.delete(self.transferred_magnetization_probs_all, i, axis=0), axis=0
         )
@@ -165,7 +166,7 @@ def _kurtosis_excluding_i(self, i: int) -> float:
             - 3
         )
 
-    def jackknife_mean(self) -> float:
+    def jackknife_mean(self) -> Real:
         r"""Compute the statistical uncertainty of the mean using the remove-one jackknife.
         If there is only one initial state (for example if $\mu = \infty$), zero uncertainty
         is returned.
@@ -175,7 +176,7 @@ def jackknife_mean(self) -> float:
         mean_i = [self._mean_excluding_i(i) for i in range(self.num_initial_states)]
         return np.std(mean_i) * np.sqrt(self.num_initial_states - 1)
 
-    def jackknife_variance(self) -> float:
+    def jackknife_variance(self) -> Real:
         r"""Compute the statistical uncertainty of the variance using the remove-one jackknife.
         If there is only one initial state (for example if $\mu = \infty$), zero uncertainty
         is returned.
@@ -187,7 +188,7 @@ def jackknife_variance(self) -> float:
         ]
         return np.std(variance_i) * np.sqrt(self.num_initial_states - 1)
 
-    def jackknife_skew(self) -> float:
+    def jackknife_skew(self) -> Real:
         r"""Compute the statistical uncertainty of the skewness using the remove-one jackknife.
         If there is only one initial state (for example if $\mu = \infty$), zero uncertainty
         is returned.
@@ -197,7 +198,7 @@ def jackknife_skew(self) -> float:
         skew_i = [self._skew_excluding_i(i) for i in range(self.num_initial_states)]
         return np.std(skew_i) * np.sqrt(self.num_initial_states - 1)
 
-    def jackknife_kurtosis(self) -> float:
+    def jackknife_kurtosis(self) -> Real:
         r"""Compute the statistical uncertainty of the kurtosis using the remove-one jackknife.
         If there is only one initial state (for example if $\mu = \infty$), zero uncertainty
         is returned.
@@ -276,35 +277,35 @@ def _transferred_magnetization(self) -> np.ndarray:
         initial = np.outer(self.num_right_initial, np.ones(num_reps, dtype=int))
         return 2 * (final - initial)
 
-    def _mean(self) -> float:
+    def _mean(self) -> Real:
         return np.mean(self.transferred_magnetization.flatten())
 
-    def _variance(self) -> float:
+    def _variance(self) -> Real:
         return np.var(self.transferred_magnetization.flatten())
 
-    def _skewness(self) -> float:
+    def _skewness(self) -> Real:
         return sstats.skew(self.transferred_magnetization.flatten())
 
-    def _kurtosis(self) -> float:
+    def _kurtosis(self) -> Real:
         return sstats.kurtosis(self.transferred_magnetization.flatten(), fisher=True)
 
-    def _mean_excluding_i(self, i: int, axis: Optional[int] = 0) -> float:
+    def _mean_excluding_i(self, i: int, axis: Optional[int] = 0) -> Real:
         tm = np.delete(self.transferred_magnetization, i, axis=axis)
         return np.mean(tm.flatten())
 
-    def _variance_excluding_i(self, i: int, axis: Optional[int] = 0) -> float:
+    def _variance_excluding_i(self, i: int, axis: Optional[int] = 0) -> Real:
         tm = np.delete(self.transferred_magnetization, i, axis=axis)
         return np.var(tm.flatten())
 
-    def _skew_excluding_i(self, i: int, axis: Optional[int] = 0) -> float:
+    def _skew_excluding_i(self, i: int, axis: Optional[int] = 0) -> Real:
         tm = np.delete(self.transferred_magnetization, i, axis=axis)
         return sstats.skew(tm.flatten())
 
-    def _kurtosis_excluding_i(self, i: int, axis: Optional[int] = 0) -> float:
+    def _kurtosis_excluding_i(self, i: int, axis: Optional[int] = 0) -> Real:
         tm = np.delete(self.transferred_magnetization, i, axis=axis)
         return sstats.kurtosis(tm.flatten(), fisher=True)
 
-    def jackknife_mean(self) -> float:
+    def jackknife_mean(self) -> Real:
         """Compute the statistical uncertainty of the mean using the remove-one jackknife.
         In the case that there is only one initial state, use the standard deviation of
         the measured transferred magnetization to estimate the uncertainty instead.
@@ -315,7 +316,7 @@ def jackknife_mean(self) -> float:
         mean_i = [self._mean_excluding_i(i) for i in range(self.num_initial_states)]
         return np.std(mean_i) * np.sqrt(self.num_initial_states - 1)
 
-    def jackknife_variance(self) -> float:
+    def jackknife_variance(self) -> Real:
         """Compute the statistical uncertainty of the variance using the remove-one jackknife.
         One initial state is removed, and the variation depending on which state is removed
         is used to estimate the uncertainty. In the case that there is only one initial state,
@@ -330,7 +331,7 @@ def jackknife_variance(self) -> float:
         variance_i = [self._variance_excluding_i(i, axis=axis) for i in range(tot)]
         return np.std(variance_i) * np.sqrt(tot - 1)
 
-    def jackknife_skew(self) -> float:
+    def jackknife_skew(self) -> Real:
         """Compute the statistical uncertainty of the skewness using the remove-one jackknife.
         One initial state is removed, and the variation depending on which state is removed
         is used to estimate the uncertainty. In the case that there is only one initial state,
@@ -345,7 +346,7 @@ def jackknife_skew(self) -> float:
         skew_i = [self._skew_excluding_i(i, axis=axis) for i in range(tot)]
         return np.std(skew_i) * np.sqrt(tot - 1)
 
-    def jackknife_kurtosis(self) -> float:
+    def jackknife_kurtosis(self) -> Real:
         """Compute the statistical uncertainty of the kurtosis using the remove-one jackknife.
         One initial state is removed, and the variation depending on which state is removed
         is used to estimate the uncertainty. In the case that there is only one initial state,
@@ -409,10 +410,10 @@ class KPZExperiment:
     def __init__(
         self,
         num_cycles: int,
-        mu: float,
+        mu: Real,
         num_init_states: int,
-        theta: float,
-        phi: float,
+        theta: Real,
+        phi: Real,
         num_qubits: Optional[Union[None, int]] = None,
     ):
         """
diff --git a/recirq/lattice_gauge/lattice_gauge_experiment.py b/recirq/lattice_gauge/lattice_gauge_experiment.py
index 5342ca58..97f5d655 100644
--- a/recirq/lattice_gauge/lattice_gauge_experiment.py
+++ b/recirq/lattice_gauge/lattice_gauge_experiment.py
@@ -13,6 +13,7 @@
 # limitations under the License.
 
 from collections.abc import Sequence
+from numbers import Real
 from typing import Any
 
 import cirq
@@ -25,7 +26,7 @@
 from recirq.lattice_gauge.lattice_gauge_grid import QubitNeighbor, LGTGrid
 
 def variational_ground_state_minimal_qubits_cols(
-    grid: LGTGrid, x_ancillary_qubits_in_cols: list[set[cirq.GridQubit]], theta: float
+    grid: LGTGrid, x_ancillary_qubits_in_cols: list[set[cirq.GridQubit]], theta: Real
 ) -> list[cirq.Moment]:
     """Moments to prepare the state from the toric code variational ansatz for two columns.
 
@@ -63,7 +64,7 @@ def variational_ground_state_minimal_qubits_cols(
 
 
 def variational_ground_state_minimal_qubits(
-    grid: LGTGrid, theta: float, extra_x_plaquette_indices: list[tuple[int, int]] = []
+    grid: LGTGrid, theta: Real, extra_x_plaquette_indices: list[tuple[int, int]] = []
 ) -> list[cirq.Moment]:
     """Moments to prepare the state from the toric code variation ansatz.
 
@@ -585,7 +586,7 @@ def x_plaquette_bitstrings(data: np.array, grid: LGTGrid) -> np.ndarray:
 
 def cnot_on_layer(
     pairs_list: Sequence[tuple[cirq.GridQubit, cirq.GridQubit]],
-    depolarization_probability: float | dict | None = None,
+    depolarization_probability: Real | dict | None = None,
 ) -> Sequence[cirq.Moment]:
     """Outputs a list of moments for CNOT between two lists, in terms of CZ gates.
 
@@ -601,7 +602,7 @@ def cnot_on_layer(
             cirq.Moment(cirq.CZ.on(qc, qt) for qc, qt in pairs_list),
             cirq.Moment(cirq.H.on_each(pair[1] for pair in pairs_list)),
         ]
-    elif isinstance(depolarization_probability, float):
+    elif isinstance(depolarization_probability, Real):
         return [
             cirq.Moment(cirq.H.on_each(pair[1] for pair in pairs_list)),
             cirq.Moment(cirq.CZ.on(qc, qt) for qc, qt in pairs_list),
diff --git a/recirq/optimize/mgd.py b/recirq/optimize/mgd.py
index 92250674..b6e2a2ac 100644
--- a/recirq/optimize/mgd.py
+++ b/recirq/optimize/mgd.py
@@ -13,6 +13,7 @@
 # limitations under the License.
 
 from typing import Callable, List, Optional, Tuple
+from numbers import Real
 
 import numpy as np
 import scipy
@@ -26,7 +27,7 @@
 
 def _get_least_squares_model_gradient(
         xs: List[np.ndarray],
-        ys: List[float],
+        ys: List[Real],
         xopt: np.ndarray,
 ) -> Tuple[np.ndarray, Pipeline]:
     """Fit a least squares quadratic model and return its gradient.
@@ -58,7 +59,7 @@ def _get_least_squares_model_gradient(
     return linear_coeffs, model
 
 
-def _random_point_in_ball(n: int, radius: float) -> np.ndarray:
+def _random_point_in_ball(n: int, radius: Real) -> np.ndarray:
     """Return a point uniformly at random from a ball centered at the origin.
 
     Args:
@@ -76,19 +77,19 @@ def _random_point_in_ball(n: int, radius: float) -> np.ndarray:
 
 
 def model_gradient_descent(
-        f: Callable[..., float],
+        f: Callable[..., Real],
         x0: np.ndarray,
         *,
         args=(),
-        rate: float = 1e-1,
-        sample_radius: float = 1e-1,
+        rate: Real = 1e-1,
+        sample_radius: Real = 1e-1,
         n_sample_points: int = 100,
-        n_sample_points_ratio: Optional[float] = None,
-        rate_decay_exponent: float = 0.0,
-        stability_constant: float = 0.0,
-        sample_radius_decay_exponent: float = 0.0,
-        tol: float = 1e-8,
-        known_values: Optional[Tuple[List[np.ndarray], List[float]]] = None,
+        n_sample_points_ratio: Optional[Real] = None,
+        rate_decay_exponent: Real = 0.0,
+        stability_constant: Real = 0.0,
+        sample_radius_decay_exponent: Real = 0.0,
+        tol: Real = 1e-8,
+        known_values: Optional[Tuple[List[np.ndarray], List[Real]]] = None,
         max_iterations: Optional[int] = None,
         max_evaluations: Optional[int] = None) -> scipy.optimize.OptimizeResult:
     """Model gradient descent algorithm for black-box optimization.
diff --git a/recirq/optimize/mgd_test.py b/recirq/optimize/mgd_test.py
index 6cd08e38..f799ba83 100644
--- a/recirq/optimize/mgd_test.py
+++ b/recirq/optimize/mgd_test.py
@@ -12,6 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+import numbers
 import numpy as np
 
 from recirq.optimize.mgd import model_gradient_descent
@@ -72,7 +73,7 @@ def test_model_gradient_descent_limited_iterations():
                                     max_iterations=15)
 
     assert isinstance(result.x, np.ndarray)
-    assert isinstance(result.fun, float)
+    assert isinstance(result.fun, numbers.Real)
     assert result.nit == 15
 
 
@@ -90,5 +91,5 @@ def test_model_gradient_descent_limited_evaluations():
                                     max_evaluations=15)
 
     assert isinstance(result.x, np.ndarray)
-    assert isinstance(result.fun, float)
+    assert isinstance(result.fun, numbers.Real)
     assert result.nfev == 12
diff --git a/recirq/optimize/minimize.py b/recirq/optimize/minimize.py
index d0f656ac..66d6c016 100644
--- a/recirq/optimize/minimize.py
+++ b/recirq/optimize/minimize.py
@@ -1,4 +1,5 @@
-from typing import Callable, Optional, Tuple
+from typing import Any, Callable, Optional, Tuple
+from numbers import Real
 
 import numpy as np
 import scipy.optimize
@@ -9,7 +10,7 @@
 OPTIMIZERS = {'mgd': model_gradient_descent, 'mpg': model_policy_gradient}
 
 
-def minimize(fun: Callable[..., float],
+def minimize(fun: Callable[..., Real],
              x0: np.ndarray,
              args: Tuple = (),
              method: Optional[str] = None,
diff --git a/recirq/optimize/mpg.py b/recirq/optimize/mpg.py
index aab3554f..efda9769 100644
--- a/recirq/optimize/mpg.py
+++ b/recirq/optimize/mpg.py
@@ -13,7 +13,8 @@
 # limitations under the License.
 
 from dataclasses import dataclass
-from typing import TYPE_CHECKING, Callable, List, Optional, Tuple
+from typing import TYPE_CHECKING, Callable, List, Optional, Tuple, Any
+from numbers import Real
 
 import numpy as np
 import scipy
@@ -31,7 +32,7 @@
 
 
 def _get_quadratic_model(
-    xs: List[np.ndarray], ys: List[float], xopt: np.ndarray
+    xs: List[np.ndarray], ys: List[Real], xopt: np.ndarray
 ) -> Pipeline:
     """Fit a least squares quadratic model.
 
@@ -61,36 +62,36 @@ def _get_quadratic_model(
 @dataclass(frozen=True)
 class _ExponentialSchedule:
     """The Exponential schedule for some hyperparameter (e.g. learning_rate)
-    
-    Exponential decay for the `learning rate`. For each `decay_steps`, the learning 
-    rate is scheduled to decay at the `decay_rate`. The `staircase` controls whether 
+
+    Exponential decay for the `learning rate`. For each `decay_steps`, the learning
+    rate is scheduled to decay at the `decay_rate`. The `staircase` controls whether
     to decay smoothly or discontinuously. After this many timesteps pass, the final
     learning rate is returned.
 
-    Args: 
-        learning rate: the initial learning rate 
-        decay_steps: the learning rate is scheduled to decay every such number of steps 
+    Args:
+        learning rate: the initial learning rate
+        decay_steps: the learning rate is scheduled to decay every such number of steps
         decay_rate: the learning rate is scheduled to decay at the such rate
-        staircase: if True, the learning rate keeps the same before every decay steps; 
-                    otherwise, the learning rate decays smoothly according 
+        staircase: if True, the learning rate keeps the same before every decay steps;
+                    otherwise, the learning rate decays smoothly according
                     to exponential interpolation.
-    
-    Returns: 
+
+    Returns:
         a class of the schedule
     """
 
-    learning_rate: float
+    learning_rate: Real
     decay_steps: int
-    decay_rate: float
+    decay_rate: Real
     staircase: bool = False
 
-    def value(self, t):
-        """Return the value of the schedule at time step t 
+    def value(self, t: Real) -> Real:
+        """Return the value of the schedule at time step t
 
-        Args: 
-            t: the time step for the schedule 
+        Args:
+            t: the time step for the schedule
 
-        Returns: 
+        Returns:
             the schedule value
         """
         m = t / self.decay_steps
@@ -106,32 +107,32 @@ def _adam_update(
     step: int,
     m: np.ndarray,
     v: np.ndarray,
-    lr_schedule=_ExponentialSchedule(0.001, 10, 0.93),
-    b1: float = 0.9,
-    b2: float = 0.999,
-    eps: float = 10 ** -8,
-):
-    """Performs a single optimization step of the optimizer Adam: a method for stochastic gradient descent  
+    lr_schedule: _ExponentialSchedule = _ExponentialSchedule(0.001, 10, 0.93),
+    b1: Real = 0.9,
+    b2: Real = 0.999,
+    eps: Real = 10 ** -8,
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Performs a single optimization step of the optimizer Adam: a method for stochastic gradient descent
 
         Adam as described in http://arxiv.org/pdf/1412.6980.pdf.
         adapted from https://github.com/HIPS/autograd/blob/master/autograd/misc/optimizers.py#L57-L71
 
         Args:
-            grad: the gradient computed at the current update.  
-            x: the current value of the parameters. 
-            step: current iteration step. 
+            grad: the gradient computed at the current update.
+            x: the current value of the parameters.
+            step: current iteration step.
             m: first moment estimate for Adam.
             v: second moment estimate for Adam.
             lr_schedule: the class of learning rate decay schedule. Defaults to ExponentialSchedule(0.001, 10, 0.93).
-            b1 (float, optional): coefficients used for computing
+            b1: coefficients used for computing
                 running averages of gradient. Defaults to 0.9.
-            b2 (float, optional): coefficients used for computing
+            b2: coefficients used for computing
                 running averages of gradient's square. Defaults to 0.999.
-            eps (float, optional): term added to the denominator to improve
+            eps: term added to the denominator to improve
                 numerical stability. Defaults to 10**-8.
 
         Returns:
-            x: the updated parameter values. 
+            x: the updated parameter values.
             m: the updated first moment estimate.
             v: the updated second moment estimate.
         """
@@ -151,28 +152,28 @@ def _adam_update(
 
 
 def model_policy_gradient(
-    f: Callable[..., float],
+    f: Callable[..., Real],
     x0: np.ndarray,
     *,
-    args=(),
-    learning_rate: float = 1e-2,
-    decay_rate: float = 0.96,
+    args: Tuple[Any, ...] = (),
+    learning_rate: Real = 1e-2,
+    decay_rate: Real = 0.96,
     decay_steps: int = 5,
-    log_sigma_init: float = -5.0,
+    log_sigma_init: Real = -5.0,
     max_iterations: int = 1000,
     batch_size: int = 10,
-    radius_coeff: float = 3.0,
+    radius_coeff: Real = 3.0,
     warmup_steps: int = 10,
     batch_size_model: int = 65536,
     save_func_vals: bool = False,
     random_state: "cirq.RANDOM_STATE_OR_SEED_LIKE" = None,
-    known_values: Optional[Tuple[List[np.ndarray], List[float]]] = None,
+    known_values: Optional[Tuple[List[np.ndarray], List[Real]]] = None,
     max_evaluations: Optional[int] = None
 ) -> scipy.optimize.OptimizeResult:
     """Model policy gradient algorithm for black-box optimization.
 
     The idea of this algorithm is to perform policy gradient, but estimate
-    the function values using a surrogate model. 
+    the function values using a surrogate model.
     The surrogate model is a least-squared quadratic
     fit to points sampled from the vicinity of the current iterate.
 
@@ -184,23 +185,23 @@ def model_policy_gradient(
         decay_rate: the learning decay rate for the Adam optimizer.
         decay_steps: the learning decay steps for the Adam optimizer.
         log_sigma_init: the initial value for the sigma of the policy
-            in the log scale. 
+            in the log scale.
         max_iterations: The maximum number of iterations to allow before
             termination.
-        batch_size: The number of points to sample in each iteration. The cost 
-            of evaluation of these samples are computed through the 
+        batch_size: The number of points to sample in each iteration. The cost
+            of evaluation of these samples are computed through the
             quantum computer cost model.
-        radius_coeff: The ratio determining the size of the radius around 
+        radius_coeff: The ratio determining the size of the radius around
             the current iterate to sample points from to build the quadratic model.
-            The ratio is with respect to the maximal ratio of the samples 
-            from the current policy. 
-        warmup_steps: The number of steps before the model policy gradient is performed. 
-            before these steps, we use the policy gradient without the model. 
-        batch_size_model: The model sample batch size. 
-            After we fit the quadratic model, we use the model to evaluate 
+            The ratio is with respect to the maximal ratio of the samples
+            from the current policy.
+        warmup_steps: The number of steps before the model policy gradient is performed.
+            before these steps, we use the policy gradient without the model.
+        batch_size_model: The model sample batch size.
+            After we fit the quadratic model, we use the model to evaluate
             on big enough batch of samples.
-        save_func_vals: whether to compute and save the function values for 
-            the current value of parameter.   
+        save_func_vals: whether to compute and save the function values for
+            the current value of parameter.
         random_state: A seed (int) or `np.random.RandomState` class to use when
             generating random values. If not set, defaults to using the module
             methods in `np.random`.
diff --git a/recirq/optimize/mpg_test.py b/recirq/optimize/mpg_test.py
index 82f41bfe..2e559b6a 100644
--- a/recirq/optimize/mpg_test.py
+++ b/recirq/optimize/mpg_test.py
@@ -13,6 +13,7 @@
 # limitations under the License.
 
 import numpy as np
+import numbers
 
 from recirq.optimize.mpg import model_policy_gradient
 
@@ -81,7 +82,7 @@ def test_model_policy_gradient_limited_iterations():
     )
 
     assert isinstance(result.x, np.ndarray)
-    assert isinstance(result.fun, float)
+    assert isinstance(result.fun, numbers.Real)
     assert result.nit == 15
 
 
@@ -103,7 +104,7 @@ def test_model_policy_gradient_limited_evaluations():
     )
 
     assert isinstance(result.x, np.ndarray)
-    assert isinstance(result.fun, float)
+    assert isinstance(result.fun, numbers.Real)
     assert result.nfev == 91
 
 
diff --git a/recirq/qcqmc/optimize_wf.py b/recirq/qcqmc/optimize_wf.py
index 47308b7f..38e2deae 100644
--- a/recirq/qcqmc/optimize_wf.py
+++ b/recirq/qcqmc/optimize_wf.py
@@ -15,6 +15,7 @@
 
 import copy
 import itertools
+import numbers
 from typing import Dict, List, Mapping, Optional, Sequence, Tuple
 
 import cirq
@@ -56,7 +57,7 @@ def get_and_check_energy(
     two_body_params: np.ndarray,
     one_body_basis_change_mat: np.ndarray,
     params: trial_wf.PerfectPairingPlusTrialWavefunctionParams,
-) -> Tuple[float, float]:
+) -> Tuple[numbers.Real, numbers.Real]:
     """Compute the energy of the ansatz circuit and check against known values where possible.
 
     Args:
@@ -225,7 +226,7 @@ def get_rotated_hamiltonians(
     one_body_basis_change_mat: np.ndarray,
     mode_qubit_map: Mapping[fermion_mode.FermionicMode, cirq.Qid],
     ordered_qubits: Sequence[cirq.Qid],
-) -> Tuple[fqe_hams.RestrictedHamiltonian, float, scipy.sparse.csc_matrix]:
+) -> Tuple[fqe_hams.RestrictedHamiltonian, numbers.Real, scipy.sparse.csc_matrix]:
     """A helper method that gets the hamiltonians in the basis of the trial_wf.
 
     Args:
@@ -267,10 +268,10 @@ def get_energy_and_check_sanity(
     unrotated_fqe_wf: fqe_wfn.Wavefunction,
     fqe_ham: fqe_hams.RestrictedHamiltonian,
     sparse_ham: scipy.sparse.csc_matrix,
-    e_core: float,
+    e_core: numbers.Real,
     mode_qubit_map: Mapping[fermion_mode.FermionicMode, cirq.Qid],
     ordered_qubits: Sequence[cirq.Qid],
-) -> float:
+) -> numbers.Real:
     """A method that checks for consistency and returns the ansatz energy.
 
     Args:
@@ -295,7 +296,7 @@ def get_energy_and_check_sanity(
     ansatz_energy = np.real_if_close(
         (np.conj(circuit_wf) @ sparse_ham @ circuit_wf)
     ).item()
-    assert isinstance(ansatz_energy, float)
+    assert isinstance(ansatz_energy, numbers.Real)
 
     fqe_energy = np.real(fqe_wf.expectationValue(fqe_ham) + e_core)
     np.testing.assert_array_almost_equal(ansatz_energy, fqe_energy)
@@ -435,8 +436,8 @@ def evaluate_energy_and_gradient(
     two_body_params: np.ndarray,
     gate_generators: List[of.FermionOperator],
     restricted: bool,
-    e_core: float,
-) -> Tuple[float, np.ndarray]:
+    e_core: numbers.Real,
+) -> Tuple[numbers.Real, np.ndarray]:
     """Evaluate gradient and cost function for optimization.
 
     Uses the linear scaling algorithm (see algo 1 from
@@ -542,7 +543,7 @@ def compute_finite_difference_grad(
     two_body_params: np.ndarray,
     ham: fqe_hams.RestrictedHamiltonian,
     initial_wf: fqe_wfn.Wavefunction,
-    dtheta: float = 1e-4,
+    dtheta: numbers.Real = 1e-4,
     restricted: bool = False,
 ):
     """Compute the parameter gradient using finite differences.
@@ -620,9 +621,9 @@ def objective(
     n_orb: int,
     restricted: bool,
     initial_orbital_rotation: np.ndarray,
-    e_core: float,
+    e_core: numbers.Real,
     do_print: bool = False,
-) -> float:
+) -> numbers.Real:
     """Helper function to compute energy from the variational parameters.
 
     Args:
@@ -671,9 +672,9 @@ def objective_and_gradient(
     gate_generators: List[of.FermionOperator],
     n_orb: int,
     restricted: bool,
-    e_core: float,
+    e_core: numbers.Real,
     do_print: bool = False,
-) -> Tuple[float, np.array]:
+) -> Tuple[numbers.Real, np.array]:
     """Helper function to compute energy and gradient from the variational parameters
 
     Args:
@@ -719,12 +720,12 @@ def optimize_parameters(
     n_one_body_params: int,
     n_two_body_params: int,
     fqe_ham: fqe_hams.RestrictedHamiltonian,
-    e_core: float,
+    e_core: numbers.Real,
     initial_orbital_rotation: Optional[np.ndarray] = None,
     restricted: bool = False,
     use_fast_gradients: bool = False,
     n_optimization_restarts: int = 1,
-    random_parameter_scale: float = 1.0,
+    random_parameter_scale: numbers.Real = 1.0,
     do_print: bool = True,
 ) -> Optional[scipy.optimize.OptimizeResult]:
     """Optimize the cost function (total energy) for the PP+ ansatz.
@@ -821,7 +822,7 @@ def get_pp_plus_params(
     *,
     hamiltonian_data: hamiltonian.HamiltonianData,
     restricted: bool = False,
-    random_parameter_scale: float = 1.0,
+    random_parameter_scale: numbers.Real = 1.0,
     initial_orbital_rotation: Optional[np.ndarray] = None,
     heuristic_layers: Tuple[layer_spec.LayerSpec, ...],
     do_pp: bool = True,
diff --git a/recirq/qml_lfe/circuit_blocks.py b/recirq/qml_lfe/circuit_blocks.py
index 5fa4be49..e36b142f 100644
--- a/recirq/qml_lfe/circuit_blocks.py
+++ b/recirq/qml_lfe/circuit_blocks.py
@@ -14,7 +14,7 @@
 
 """Sycamore native circuit building blocks."""
 from typing import Callable, Dict, List, Tuple
-from numbers import Number
+from numbers import Real
 
 import os
 import cirq
@@ -29,7 +29,7 @@ def _load_circuit(fname: str) -> cirq.Circuit:
         return cirq.read_json(f)
 
 
-def tsym_block(qubits: List[cirq.Qid], params: List[Number]) -> List[cirq.Operation]:
+def tsym_block(qubits: List[cirq.Qid], params: List[Real]) -> List[cirq.Operation]:
     """Create a tsym block.
 
     Y axis rotations followed by a real entried permutation operation.
@@ -48,7 +48,7 @@ def tsym_block(qubits: List[cirq.Qid], params: List[Number]) -> List[cirq.Operat
     return rots + list(mapped_circuit.all_operations())
 
 
-def _get_op(code: Number) -> cirq.Gate:
+def _get_op(code: Real) -> cirq.Gate:
     if code <= 4 / 3:
         return cirq.X ** 0.5
     elif code <= 2 * (4 / 3):
@@ -60,7 +60,7 @@ def _get_op(code: Number) -> cirq.Gate:
 
 
 def scrambling_block(
-    qubits: List[cirq.Qid], params: List[Number]
+    qubits: List[cirq.Qid], params: List[Real]
 ) -> List[cirq.Operation]:
     """Create a scrambling block.
 
@@ -84,7 +84,7 @@ def block_1d_circuit(
     qubits: List[cirq.Qid],
     depth: int,
     block_fn: Callable,
-    random_source: List[List[Number]],
+    random_source: List[List[Real]],
 ) -> cirq.Circuit:
     """Create a 1D block structure circuit using block_fn.