diff --git a/.github/workflows/regression-tests.yml b/.github/workflows/regression-tests.yml
index 7b960e9..5f5af7e 100644
--- a/.github/workflows/regression-tests.yml
+++ b/.github/workflows/regression-tests.yml
@@ -16,7 +16,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
 
     steps:
     - uses: actions/checkout@v3
@@ -24,20 +24,19 @@ jobs:
       uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python-version }}
-    - name: Install dependencies
+    - name: Install base dependencies
       run: |
         python -m pip install --upgrade pip
+        pip install .
+        # Basic check for minimal deployed env requirements
+        python -c "import pyttb"
+    - name: Install dev dependencies
+      run: |
         python -m pip install --upgrade coverage coveralls sphinx_rtd_theme
         pip install ".[dev]"
     - name: Check auto-formatters
       run: |
         black --check .
-#    - name: Lint with flake8
-#      run: |
-#        # stop the build if there are Python syntax errors or undefined names
-#        flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
-#        # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
-#        flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Run tests
       run: |
         coverage run --source pyttb -m pytest tests/
diff --git a/pyttb/gcp/handles.py b/pyttb/gcp/handles.py
index c95ceab..06a8dff 100644
--- a/pyttb/gcp/handles.py
+++ b/pyttb/gcp/handles.py
@@ -115,15 +115,17 @@ def huber_grad(data: ttb.tensor, model: ttb.tensor, threshold: float) -> np.ndar
     ) * np.logical_not(below_threshold)
 
 
+# FIXME: Num trials should be enforced as integer here and in MATLAB
+# requires updating our regression test values to calculate MATLAB integer version
 def negative_binomial(
-    data: np.ndarray, model: np.ndarray, num_trials: int
+    data: np.ndarray, model: np.ndarray, num_trials: float
 ) -> np.ndarray:
     """Return objective function for negative binomial distributions"""
     return (num_trials + data) * np.log(model + 1) - data * np.log(model + EPS)
 
 
 def negative_binomial_grad(
-    data: np.ndarray, model: np.ndarray, num_trials: int
+    data: np.ndarray, model: np.ndarray, num_trials: float
 ) -> np.ndarray:
     """Return gradient function for negative binomial distributions"""
     return (num_trials + 1) / (1 + model) - data / (model + EPS)
diff --git a/pyttb/gcp/samplers.py b/pyttb/gcp/samplers.py
index 161571a..d7dbda0 100644
--- a/pyttb/gcp/samplers.py
+++ b/pyttb/gcp/samplers.py
@@ -349,7 +349,7 @@ def zeros(
 
     # Select out just the zeros
     tmpidx = tt_sub2ind(data.shape, tmpsubs)
-    iszero = np.logical_not(np.in1d(tmpidx, nz_idx))
+    iszero = np.logical_not(np.isin(tmpidx, nz_idx))
     tmpsubs = tmpsubs[iszero, :]
 
     # Trim back to desired numb of samples
@@ -425,7 +425,7 @@ def semistrat(data: ttb.sptensor, num_nonzeros: int, num_zeros: int) -> sample_t
 
 
 def stratified(
-    data: ttb.sptensor,
+    data: Union[ttb.sptensor, ttb.tensor],
     nz_idx: np.ndarray,
     num_nonzeros: int,
     num_zeros: int,
@@ -450,6 +450,9 @@ def stratified(
     -------
     Subscripts, values, and weights of samples (Nonzeros then zeros).
     """
+    assert isinstance(
+        data, ttb.sptensor
+    ), "For stratified sampling Sparse Tensor must be provided"
     [nonzero_subs, nonzero_vals] = nonzeros(data, num_nonzeros, with_replacement=True)
     nonzero_weights = np.ones((num_nonzeros,))
     if num_nonzeros > 0:
diff --git a/pyttb/ktensor.py b/pyttb/ktensor.py
index 475f892..38a0a25 100644
--- a/pyttb/ktensor.py
+++ b/pyttb/ktensor.py
@@ -22,13 +22,13 @@
 
 import numpy as np
 import scipy.sparse.linalg
-from typing_extensions import Self
 
 import pyttb as ttb
 from pyttb.pyttb_utils import (
     get_mttkrp_factors,
     isrow,
     isvector,
+    np_to_python,
     tt_dimscheck,
     tt_ind2sub,
 )
@@ -695,7 +695,8 @@ def extract(
                         invalid_entries.append(component)
                 if len(invalid_entries) > 0:
                     assert False, (
-                        f"Invalid component indices to be extracted: {invalid_entries} "
+                        f"Invalid component indices to be extracted: "
+                        f"{np_to_python(invalid_entries)} "
                         f"not in range({self.ncomponents})"
                     )
                 new_weights = self.weights[components]
@@ -706,7 +707,7 @@ def extract(
         else:
             assert False, "Input parameter must be an int, tuple, list or numpy.ndarray"
 
-    def fixsigns(self, other: Optional[ktensor] = None) -> Self:  # noqa: PLR0912
+    def fixsigns(self, other: Optional[ktensor] = None) -> ktensor:  # noqa: PLR0912
         """
         Change the elements of a :class:`pyttb.ktensor` in place so that the
         largest magnitude entries for each column vector in each factor
@@ -1184,7 +1185,7 @@ def normalize(
         sort: Optional[bool] = False,
         normtype: float = 2,
         mode: Optional[int] = None,
-    ) -> Self:
+    ) -> ktensor:
         """
         Normalize the columns of the factor matrices of a
         :class:`pyttb.ktensor` in place, then optionally
@@ -1405,7 +1406,7 @@ def permute(self, order: np.ndarray) -> ktensor:
 
         return ttb.ktensor([self.factor_matrices[i] for i in order], self.weights)
 
-    def redistribute(self, mode: int) -> Self:
+    def redistribute(self, mode: int) -> ktensor:
         """
         Distribute weights of a :class:`pyttb.ktensor` to the specified mode.
         The redistribution is performed in place.
@@ -1621,7 +1622,7 @@ def score(
 
             # Rearrange the components of A according to the best matching
             foo = np.arange(RA)
-            tf = np.in1d(foo, best_perm)
+            tf = np.isin(foo, best_perm)
             best_perm[RB : RA + 1] = foo[~tf]
             A.arrange(permutation=best_perm)
             return best_score, A, flag, best_perm
@@ -1999,7 +2000,7 @@ def ttv(
             factor_matrices.append(self.factor_matrices[i])
         return ttb.ktensor(factor_matrices, new_weights, copy=False)
 
-    def update(self, modes: Union[int, Iterable[int]], data: np.ndarray) -> Self:
+    def update(self, modes: Union[int, Iterable[int]], data: np.ndarray) -> ktensor:
         """
         Updates a :class:`pyttb.ktensor` in the specific dimensions with the
         values in `data` (in vector or matrix form). The value of `modes` must
diff --git a/pyttb/pyttb_utils.py b/pyttb/pyttb_utils.py
index 2a467b6..84bace2 100644
--- a/pyttb/pyttb_utils.py
+++ b/pyttb/pyttb_utils.py
@@ -7,7 +7,16 @@
 
 from enum import Enum
 from inspect import signature
-from typing import List, Literal, Optional, Tuple, Union, get_args, overload
+from typing import (
+    Iterable,
+    List,
+    Literal,
+    Optional,
+    Tuple,
+    Union,
+    get_args,
+    overload,
+)
 
 import numpy as np
 
@@ -921,3 +930,22 @@ def gather_wrap_dims(
 
     assert rdims is not None and cdims is not None
     return rdims.astype(int), cdims.astype(int)
+
+
+def np_to_python(
+    iterable: Iterable,
+) -> Iterable:
+    """Convert a structure containing numpy scalars to pure python types.
+
+    Mostly useful for prettier printing post numpy 2.0.
+
+    Parameters
+    ----------
+    iterable:
+        Structure potentially containing numpy scalars.
+    """
+    output_type = type(iterable)
+    return output_type(  # type: ignore [call-arg]
+        element.item() if isinstance(element, np.generic) else element
+        for element in iterable
+    )
diff --git a/pyttb/sptenmat.py b/pyttb/sptenmat.py
index 8e8f647..5f99b39 100644
--- a/pyttb/sptenmat.py
+++ b/pyttb/sptenmat.py
@@ -129,7 +129,7 @@ def __init__(  # noqa: PLR0913
             # Sum the corresponding values
             # Squeeze to convert from column vector to row vector
             newvals = accumarray(
-                loc, np.squeeze(vals, axis=1), size=newsubs.shape[0], func=sum
+                loc.flatten(), np.squeeze(vals, axis=1), size=newsubs.shape[0], func=sum
             )
 
         # Find the nonzero indices of the new values
diff --git a/pyttb/sptensor.py b/pyttb/sptensor.py
index de1f4d2..28690a7 100644
--- a/pyttb/sptensor.py
+++ b/pyttb/sptensor.py
@@ -33,6 +33,7 @@
     gather_wrap_dims,
     get_index_variant,
     get_mttkrp_factors,
+    np_to_python,
     tt_dimscheck,
     tt_ind2sub,
     tt_intersect_rows,
@@ -325,7 +326,10 @@ def from_aggregator(
             # Sum the corresponding values
             # Squeeze to convert from column vector to row vector
             newvals = accumarray(
-                loc, np.squeeze(vals), size=newsubs.shape[0], func=function_handle
+                loc.flatten(),
+                np.squeeze(vals),
+                size=newsubs.shape[0],
+                func=function_handle,
             )
 
         # Find the nonzero indices of the new values
@@ -445,7 +449,10 @@ def collapse(
 
         # Check for the case where we accumulate over *all* dimensions
         if remdims.size == 0:
-            return function_handle(self.vals.transpose()[0])
+            result = function_handle(self.vals.transpose()[0])
+            if isinstance(result, np.generic):
+                result = result.item()
+            return result
 
         # Calculate the size of the result
         newsize = np.array(self.shape)[remdims]
@@ -1319,7 +1326,7 @@ def nnz(self) -> int:
             return 0
         return self.subs.shape[0]
 
-    def norm(self) -> np.floating:
+    def norm(self) -> float:
         """
         Compute the norm (i.e., Frobenius norm, or square root of the sum of
         squares of entries) of the :class:`pyttb.sptensor`.
@@ -1339,7 +1346,7 @@ def norm(self) -> np.floating:
         >>> S.norm() # doctest: +ELLIPSIS
         5.47722557...
         """
-        return np.linalg.norm(self.vals)
+        return np.linalg.norm(self.vals).item()
 
     def nvecs(self, n: int, r: int, flipsign: bool = True) -> np.ndarray:
         """
@@ -1945,7 +1952,7 @@ def ttv(  # noqa: PLR0912
 
         # Case 0: If all dimensions were used, then just return the sum
         if remdims.size == 0:
-            return np.sum(newvals)
+            return np.sum(newvals).item()
 
         # Otherwise, figure out new subscripts and accumulate the results.
         newsiz = np.array(self.shape, dtype=int)[remdims]
@@ -3439,10 +3446,10 @@ def __repr__(self):  # pragma: no cover
         """
         nz = self.nnz
         if nz == 0:
-            s = f"empty sparse tensor of shape {self.shape!r}"
+            s = f"empty sparse tensor of shape {np_to_python(self.shape)!r}"
             return s
 
-        s = f"sparse tensor of shape {self.shape!r}"
+        s = f"sparse tensor of shape {np_to_python(self.shape)!r}"
         s += f" with {nz} nonzeros\n"
 
         # Stop insane printouts
diff --git a/pyttb/sumtensor.py b/pyttb/sumtensor.py
index c0683d5..b1cbc1a 100644
--- a/pyttb/sumtensor.py
+++ b/pyttb/sumtensor.py
@@ -13,6 +13,7 @@
 import numpy as np
 
 import pyttb as ttb
+from pyttb.pyttb_utils import np_to_python
 
 
 class sumtensor:
@@ -115,7 +116,10 @@ def __repr__(self):
         """
         if len(self.parts) == 0:
             return "Empty sumtensor"
-        s = f"sumtensor of shape {self.shape} with {len(self.parts)} parts:"
+        s = (
+            f"sumtensor of shape {np_to_python(self.shape)} "
+            f"with {len(self.parts)} parts:"
+        )
         for i, part in enumerate(self.parts):
             s += f"\nPart {i}: \n"
             s += indent(str(part), prefix="\t")
@@ -298,15 +302,15 @@ def innerprod(
 
         Examples
         --------
-        >>> T1 = ttb.tensor(np.array([[1, 0], [0, 4]]))
+        >>> T1 = ttb.tensor(np.array([[1., 0.], [0., 4.]]))
         >>> T2 = T1.to_sptensor()
         >>> S = ttb.sumtensor([T1, T2])
         >>> T1.innerprod(T1)
-        17
+        17.0
         >>> T1.innerprod(T2)
-        17
+        17.0
         >>> S.innerprod(T1)
-        34
+        34.0
         """
         result = self.parts[0].innerprod(other)
         for part in self.parts[1:]:
diff --git a/pyttb/tenmat.py b/pyttb/tenmat.py
index 24d3cab..64a2eea 100644
--- a/pyttb/tenmat.py
+++ b/pyttb/tenmat.py
@@ -30,7 +30,7 @@ def __init__(  # noqa: PLR0912
     ):
         """
         Construct a :class:`pyttb.tenmat` from explicit components.
-        If you already have a tensor see :method:`pyttb.tensor.to_tenmat`.
+        If you already have a tensor see :meth:`pyttb.tensor.to_tenmat`.
 
         Parameters
         ----------
@@ -176,9 +176,12 @@ def copy(self) -> tenmat:
         >>> TM2 = TM1
         >>> TM3 = TM1.copy()
         >>> TM1[0,0] = 3
-        >>> TM1[0,0] == TM2[0,0]
+
+        # Item to convert numpy boolean to python boolena for nicer printing
+
+        >>> (TM1[0,0] == TM2[0,0]).item()
         True
-        >>> TM1[0,0] == TM3[0,0]
+        >>> (TM1[0,0] == TM3[0,0]).item()
         False
         """
         # Create tenmat
diff --git a/pyttb/tensor.py b/pyttb/tensor.py
index b02b40e..c3da3e8 100644
--- a/pyttb/tensor.py
+++ b/pyttb/tensor.py
@@ -23,6 +23,7 @@
     gather_wrap_dims,
     get_index_variant,
     get_mttkrp_factors,
+    np_to_python,
     tt_dimscheck,
     tt_ind2sub,
     tt_sub2ind,
@@ -254,7 +255,10 @@ def collapse(
 
         # Check for the case where we accumulate over *all* dimensions
         if remdims.size == 0:
-            return fun(self.data.flatten("F"))
+            result = fun(self.data.flatten("F"))
+            if isinstance(result, np.generic):
+                result = result.item()
+            return result
 
         ## Calculate the shape of the result
         newshape = tuple(np.array(self.shape)[remdims])
@@ -320,7 +324,7 @@ def contract(self, i1: int, i2: int) -> Union[np.ndarray, tensor]:
 
         # Easy case - returns a scalar
         if self.ndims == 2:
-            return np.trace(self.data)
+            return np.trace(self.data).item()
 
         # Remaining dimensions after trace
         remdims = np.setdiff1d(np.arange(0, self.ndims), np.array([i1, i2])).astype(int)
@@ -546,9 +550,9 @@ def to_tenmat(
         # Verify inputs
         if rdims is None and cdims is None:
             assert False, "Either rdims or cdims or both must be specified."
-        if rdims is not None and not sum(np.in1d(rdims, alldims)) == len(rdims):
+        if rdims is not None and not sum(np.isin(rdims, alldims)) == len(rdims):
             assert False, "Values in rdims must be in [0, source.ndims]."
-        if cdims is not None and not sum(np.in1d(cdims, alldims)) == len(cdims):
+        if cdims is not None and not sum(np.isin(cdims, alldims)) == len(cdims):
             assert False, "Values in cdims must be in [0, source.ndims]."
 
         rdims, cdims = gather_wrap_dims(n, rdims, cdims, cdims_cyclic)
@@ -587,19 +591,19 @@ def innerprod(
 
         Examples
         --------
-        >>> T = ttb.tensor(np.array([[1, 0], [0, 4]]))
+        >>> T = ttb.tensor(np.array([[1., 0.], [0., 4.]]))
         >>> T.innerprod(T)
-        17
+        17.0
         >>> S = T.to_sptensor()
         >>> T.innerprod(S)
-        17
+        17.0
         """
         if isinstance(other, ttb.tensor):
             if self.shape != other.shape:
                 assert False, "Inner product must be between tensors of the same size"
             x = np.reshape(self.data, (self.data.size,), order="F")
             y = np.reshape(other.data, (other.data.size,), order="F")
-            return x.dot(y)
+            return x.dot(y).item()
         if isinstance(other, (ttb.ktensor, ttb.sptensor, ttb.ttensor)):
             # Reverse arguments and call specializer code
             return other.innerprod(self)
@@ -1018,7 +1022,7 @@ def norm(self) -> float:
         # default of np.linalg.norm is to vectorize the data and compute the vector
         # norm, which is equivalent to the Frobenius norm for multidimensional arrays.
         # However, the argument 'fro' only works for 1-D and 2-D arrays currently.
-        return float(np.linalg.norm(self.data))
+        return np.linalg.norm(self.data).item()
 
     def nvecs(self, n: int, r: int, flipsign: bool = True) -> np.ndarray:
         """
@@ -1670,7 +1674,7 @@ def ttv(
         # If needed, convert the final result back to tensor
         if n > 0:
             return ttb.tensor(c, tuple(sz[0:n]), copy=False)
-        return c[0]
+        return c[0].item()
 
     def ttsv(
         self,
@@ -1978,7 +1982,7 @@ def __getitem__(self, item):  # noqa: PLR0912
             # If the size is zero, then the result is returned as a scalar
             # otherwise, we convert the result to a tensor
             if newsiz.size == 0:
-                a = newdata
+                a = newdata.item()
             else:
                 a = ttb.tensor(newdata, copy=False)
             return a
@@ -2508,7 +2512,7 @@ def __repr__(self):
             return s
 
         s = ""
-        s += f"tensor of shape {self.shape}"
+        s += f"tensor of shape {np_to_python(self.shape)}"
 
         if self.ndims == 1:
             s += "\ndata"
diff --git a/pyttb/ttensor.py b/pyttb/ttensor.py
index 8767f02..7a4d42f 100644
--- a/pyttb/ttensor.py
+++ b/pyttb/ttensor.py
@@ -106,9 +106,12 @@ def copy(self) -> ttensor:
         >>> second = first
         >>> third = second.copy()
         >>> first.factor_matrices[0][0,0] = 2
-        >>> first.factor_matrices[0][0,0] == second.factor_matrices[0][0,0]
+
+        # Item to convert numpy boolean to python boolena for nicer printing
+
+        >>> (first.factor_matrices[0][0,0] == second.factor_matrices[0][0,0]).item()
         True
-        >>> first.factor_matrices[0][0,0] == third.factor_matrices[0][0,0]
+        >>> (first.factor_matrices[0][0,0] == third.factor_matrices[0][0,0]).item()
         False
         """
         return ttb.ttensor(self.core, self.factor_matrices, copy=True)
diff --git a/tests/test_tenmat.py b/tests/test_tenmat.py
index ad8db71..04a426f 100644
--- a/tests/test_tenmat.py
+++ b/tests/test_tenmat.py
@@ -355,7 +355,7 @@ def test_tenmat_ctranspose(sample_tenmat_4way):
 def test_tenmat_double(sample_tenmat_4way):
     (params, tenmatInstance) = sample_tenmat_4way
 
-    assert (tenmatInstance.double() == tenmatInstance.data.astype(np.float_)).all()
+    assert (tenmatInstance.double() == tenmatInstance.data.astype(np.float64)).all()
 
 
 def test_tenmat_ndims(sample_tenmat_4way):