diff --git a/.github/workflows/actions.yml b/.github/workflows/actions.yml
index 6a10cdb1..dd0fab1e 100644
--- a/.github/workflows/actions.yml
+++ b/.github/workflows/actions.yml
@@ -9,7 +9,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.9, '3.10', '3.11', '3.12']
+        python-version: ['3.11', '3.12', '3.13', '3.14']
 
     steps:
     - uses: actions/checkout@v4
diff --git a/meson.build b/meson.build
index e0e877fa..26c18f60 100644
--- a/meson.build
+++ b/meson.build
@@ -1,7 +1,7 @@
 project(
     'scikit-digital-health',
     'c',
-    version: '0.17.9',
+    version: '0.17.10',
     license: 'MIT',
     meson_version: '>=1.1',
 )
diff --git a/pyproject.toml b/pyproject.toml
index 5c2016be..b1686565 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -47,9 +47,9 @@ classifiers = [
 
 requires-python = ">=3.9"
 dependencies = [
-    "numpy>=1.25.0",
+    "numpy>=2.0.0",
     "scipy>=1.12.0",
-    "pandas>=1.0.0",
+    "pandas>=2.0.0",
     "lightgbm>=2.3.0",
     "pywavelets",
     "scikit-learn",
diff --git a/src/skdh/activity/core.py b/src/skdh/activity/core.py
index 112036fa..7b1d4bb6 100644
--- a/src/skdh/activity/core.py
+++ b/src/skdh/activity/core.py
@@ -27,7 +27,7 @@
     isnan,
 )
 from numpy.linalg import norm
-from pandas import Timedelta, Timestamp
+from pandas import Timedelta, Timestamp, DataFrame
 import matplotlib
 from matplotlib.backends.backend_pdf import PdfPages
 import matplotlib.lines as mlines
@@ -82,6 +82,8 @@ class ActivityLevelClassification(BaseProcess):
         Two (2) element array-like of the base and period of the window to use for
         determining days. Default is (0, 24), which will look for days starting at
         midnight and lasting 24 hours. None removes any day-based windowing.
+    save_epoch_data : bool
+        Save the epochs with calculated metrics. Default is False.
 
     Notes
     -----
@@ -137,6 +139,7 @@ def __init__(
         min_wear_time=10,
         cutpoints="migueles_wrist_adult",
         day_window=(0, 24),
+        save_epoch_data=False,
     ):
         # make sure that the short_wlen is a factor of 60, and if not send it to
         # nearest factor
@@ -165,6 +168,7 @@ def __init__(
             min_wear_time=min_wear_time,
             cutpoints=cutpoints_,
             day_window=day_window,
+            save_epoch_data=save_epoch_data,
         )
 
         self.wlen = short_wlen
@@ -181,6 +185,13 @@ def __init__(
         else:
             self.day_key = tuple(day_window)
 
+        self.save_epochs = save_epoch_data
+        self.epoch_data = {
+            "time": [],
+            "metric": [],
+            "intensity": []
+        }
+
         # enable plotting as a public method
         self.setup_plotting = self._setup_plotting
         self._update_buttons = []
@@ -228,6 +239,39 @@ def __init__(
             for lvl in self.act_levels
         ]
 
+    def save_results(self, results, file_name):
+        """
+        Save the results of the processing pipeline to a csv file. Will also
+        save per epoch metric value if `save_epoch_data` is true. The file name 
+        for per epoch results is the same as the activity endpoints file with 
+        "_per_epoch_predictions" added to the end.
+
+        Parameters
+        ----------
+        results : dict
+            Dictionary of results from the output of predict
+        file_name : str
+            File name. Can be optionally formatted (see Notes)
+
+        Notes
+        -----
+        Available format variables available:
+
+        - date: todays date expressed in yyyymmdd format.
+        - name: process name.
+        - file: file name used in the pipeline, or "" if not found.
+        """
+        file_name = super().save_results(results, file_name)
+
+        if self.save_epochs:
+            file_name = Path(file_name)
+
+            new_name = file_name.stem + "_per_epoch_predictions" + file_name.suffix
+            epoch_file = file_name.with_name(new_name)
+
+            df = DataFrame(self.epoch_data)
+            df.to_csv(epoch_file, index=False)
+
     def add(self, endpoint):
         """
         Add an endpoint to the list to be calculated.
@@ -478,7 +522,7 @@ def predict(self, *, time, accel, fs=None, wear=None, tz_name=None, **kwargs):
 
             # compute waking hours activity endpoints
             self._compute_awake_activity_endpoints(
-                res, accel, fs, iday, dwear_starts, dwear_stops, nwlen, nwlen_60, epm
+                res, time, accel, fs, iday, dwear_starts, dwear_stops, nwlen, nwlen_60, epm
             )
             # compute sleeping hours activity endpoints
             self._compute_sleep_activity_endpoints(
@@ -517,7 +561,7 @@ def _initialize_awake_values(self, results, day_n):
                 results[endpt.name][day_n] = 0.0
 
     def _compute_awake_activity_endpoints(
-        self, results, accel, fs, day_n, starts, stops, n_wlen, n_wlen_60, epm
+        self, results, time, accel, fs, day_n, starts, stops, n_wlen, n_wlen_60, epm
     ):
         # initialize values from nan to 0.0. Do this here because days with less than
         # minimum hours should have nan values
@@ -538,6 +582,14 @@ def _compute_awake_activity_endpoints(
                 accel[start:stop], n_wlen_60, fs, **self.cutpoints["kwargs"]
             )
 
+            # handle saving the epoch data
+            self._handle_epoch_data(
+                n_wlen,
+                fs,
+                time[start:stop],
+                acc_metric,
+            )
+
             for endpoint in self.wake_endpoints:
                 endpoint.predict(
                     results, day_n, acc_metric, acc_metric_60, self.wlen, epm
@@ -740,6 +792,30 @@ def _finalize_plots(self):  # pragma: no cover
 
         pp.close()
 
+    def _handle_epoch_data(self, wlen, fs, time, metric):
+        if self.save_epochs:
+            self.epoch_data["time"].extend(
+                time[:-wlen+1:wlen]
+            )
+            self.epoch_data["metric"].extend(
+                metric
+            )
+            # handle the threshold/intensity stuff
+            intensity = full(metric.size, "sedentary")
+
+            light_lt, light_ut = get_level_thresholds("light", self.cutpoints)
+            mod_lt, mod_ut = get_level_thresholds("moderate", self.cutpoints)
+            vig_lt, _ = get_level_thresholds("vig", self.cutpoints)
+
+            light_mask = (metric >= light_lt) & (metric < light_ut)
+            mod_mask = (metric >= mod_lt) & (metric < mod_ut)
+            vig_mask = metric >= vig_lt
+
+            intensity[light_mask] = "light"
+            intensity[mod_mask] = "moderate"
+            intensity[vig_mask] = "vigorous"
+
+            self.epoch_data['intensity'].extend(intensity)
 
 class StaudenmayerClassification(BaseProcess):
     """
diff --git a/src/skdh/features/core.py b/src/skdh/features/core.py
index 72c8a15a..512690fb 100644
--- a/src/skdh/features/core.py
+++ b/src/skdh/features/core.py
@@ -262,7 +262,7 @@ def compute(
         # standardize the input signal
         if isinstance(signal, DataFrame):
             columns = columns if columns is not None else signal.columns
-            x = signal[columns].values.astype(float64)
+            x = signal[columns].to_numpy(copy=True).astype(float64)
         else:
             try:
                 x = asarray(signal, dtype=float64)
diff --git a/src/skdh/io/base.py b/src/skdh/io/base.py
index 2f542162..129db6c3 100644
--- a/src/skdh/io/base.py
+++ b/src/skdh/io/base.py
@@ -106,7 +106,7 @@ def handle_naive_timestamps(time, is_local, tz_name=None):
                 .total_seconds()
             )
 
-            time += offset
+            time = time + offset
         else:  # is_local, tz_name is None
             warn(
                 "Timestamps are local but naive, and no time-zone information is available. "
diff --git a/src/skdh/io/csv.py b/src/skdh/io/csv.py
index cb6ca2d4..f09b26f0 100644
--- a/src/skdh/io/csv.py
+++ b/src/skdh/io/csv.py
@@ -20,6 +20,7 @@
     allclose,
     nonzero,
     diff,
+    int64,
 )
 from pandas import read_csv, to_datetime
 
@@ -266,7 +267,7 @@ def handle_timestamp_inconsistency_np(self, fill_dict, time, fs, data):
             t_delta = tile(arange(0, 1, 1 / n_samples), int(n_blocks))
 
             # add the time delta so that we have unique timestamps
-            time += t_delta
+            time = time + t_delta
 
         # check if we are filling gaps or not
         if self.fill_gaps:
@@ -332,29 +333,29 @@ def predict(self, *, file, tz_name=None, **kwargs):
         # update the to_datetime_kwargs based on tz_name.  tz_name==None (utc=False)
         self.to_datetime_kw.update({"utc": tz_name is not None})
 
-        # convert time column to a datetime column. Give a unique name so we shouldnt overwrite
-        raw[self.time_col_name] = to_datetime(
+        # convert time column to a datetime column.
+        time_series = to_datetime(
             raw[self.time_col_name], **self.to_datetime_kw
         )
 
         # convert timestamps if necessary
         if tz_name is not None:
             # convert, and then remove the timezone so its naive again, but now in local time
-            raw[self.time_col_name] = raw[self.time_col_name].dt.tz_convert(tz_name)
+            time_series = time_series.dt.tz_convert(tz_name)
+    
 
         # now handle data gaps and second level timestamps, etc
         # raw, fs = self.handle_timestamp_inconsistency(raw, fill_values)
 
         # get the time values and convert to seconds
-        time = (
-            raw[self.time_col_name].astype(int).values / 1e9
-        )  # int gives ns, convert to s
+        time = time_series.dt.as_unit('s').astype(int64).to_numpy(copy=True)
+        # first convert to 's' representation, then to int gives correct values
 
         data = {}
         # grab the data we expect
         for dstream in self.column_names:
             try:
-                data[dstream] = raw[self.column_names[dstream]].values
+                data[dstream] = raw[self.column_names[dstream]].to_numpy(copy=True)
             except KeyError:
                 warn(
                     f"Data stream {dstream} specified in column names but all "
diff --git a/src/skdh/sleep/sleep.py b/src/skdh/sleep/sleep.py
index 84f26880..e9f2d9b9 100644
--- a/src/skdh/sleep/sleep.py
+++ b/src/skdh/sleep/sleep.py
@@ -47,7 +47,7 @@ class Sleep(BaseProcess):
     min_rest_block : int, optional
         Number of minutes required to consider a rest period valid. Default is 30 minutes.
     max_activity_break : int, optional
-        Number of minutes of activity allowed to interrupt the major rest period. Default is 30
+        Number of minutes of activity allowed to interrupt the major rest period. Default is 60
         minutes.
     tso_min_thresh : float, optional
         Minimum allowed z-angle threshold for determining major rest period. Default is 0.1.
diff --git a/tests/context/conftest.py b/tests/context/conftest.py
index 458437f9..4e7f831e 100644
--- a/tests/context/conftest.py
+++ b/tests/context/conftest.py
@@ -7,7 +7,7 @@
 def ambulation_positive_data(path_tests):
     accel = read_csv(
         path_tests / "context" / "data" / "ambulation_test1.csv", header=None
-    ).values
+    ).to_numpy(copy=True)
     time = arange(0, len(accel) / 20, 1 / 20)
 
     return time, accel
@@ -17,7 +17,7 @@ def ambulation_positive_data(path_tests):
 def ambulation_negative_data(path_tests):
     accel = read_csv(
         path_tests / "context" / "data" / "ambulation_test2.csv", header=None
-    ).values
+    ).to_numpy(copy=True)
     time = arange(0, len(accel) / 20, 1 / 20)
 
     return time, accel
@@ -27,7 +27,7 @@ def ambulation_negative_data(path_tests):
 def ambulation_negative_data_50hz(path_tests):
     accel = read_csv(
         path_tests / "context" / "data" / "ambulation_test2.csv", header=None
-    ).values
+    ).to_numpy(copy=True)
     time = arange(0, len(accel) / 50, 1 / 50)
 
     return time, accel
@@ -38,7 +38,7 @@ def motion_positive_data_100hz(path_tests):
     fs = 100
     accel = read_csv(
         path_tests / "context" / "data" / "motion_test_positive.csv"
-    ).values
+    ).to_numpy(copy=True)
     accel = vstack([accel, accel, accel])
     time = arange(0, len(accel) / fs, 1 / fs)
 
@@ -50,7 +50,7 @@ def motion_negative_data_100hz(path_tests):
     fs = 100
     accel = read_csv(
         path_tests / "context" / "data" / "motion_test_negative.csv"
-    ).values
+    ).to_numpy(copy=True)
     time = arange(0, len(accel) / fs, 1 / fs)
 
     return time, accel, fs
@@ -61,7 +61,7 @@ def motion_positive_data_20hz(path_tests):
     fs = 20
     accel = read_csv(
         path_tests / "context" / "data" / "motion_test_positive.csv"
-    ).values
+    ).to_numpy(copy=True)
     accel = vstack([accel, accel, accel])
     time = arange(0, len(accel) / fs, 1 / fs)
 
@@ -73,7 +73,7 @@ def motion_negative_data_20hz(path_tests):
     fs = 20
     accel = read_csv(
         path_tests / "context" / "data" / "motion_test_negative.csv"
-    ).values
+    ).to_numpy(copy=True)
     time = arange(0, len(accel) / fs, 1 / fs)
 
     return time, accel, fs
diff --git a/tests/gait/data/gait_results2_apcwt.npz b/tests/gait/data/gait_results2_apcwt.npz
index f4c962ee..e7d0ff04 100644
Binary files a/tests/gait/data/gait_results2_apcwt.npz and b/tests/gait/data/gait_results2_apcwt.npz differ
diff --git a/tests/gait/data/gait_results2_vcwt.npz b/tests/gait/data/gait_results2_vcwt.npz
index 51d90f87..8cf703d1 100644
Binary files a/tests/gait/data/gait_results2_vcwt.npz and b/tests/gait/data/gait_results2_vcwt.npz differ
diff --git a/tests/gait/data/gait_results_apcwt.npz b/tests/gait/data/gait_results_apcwt.npz
index aaa8c4b7..24caeabe 100644
Binary files a/tests/gait/data/gait_results_apcwt.npz and b/tests/gait/data/gait_results_apcwt.npz differ
diff --git a/tests/gait/data/gait_results_vcwt.npz b/tests/gait/data/gait_results_vcwt.npz
index e571f032..775612ba 100644
Binary files a/tests/gait/data/gait_results_vcwt.npz and b/tests/gait/data/gait_results_vcwt.npz differ
diff --git a/tests/gait/data/manual_gait_results_apcwt.xlsx b/tests/gait/data/manual_gait_results_apcwt.xlsx
index 4bb3eb90..2ab9291b 100644
Binary files a/tests/gait/data/manual_gait_results_apcwt.xlsx and b/tests/gait/data/manual_gait_results_apcwt.xlsx differ
diff --git a/tests/gait/data/manual_gait_results_vcwt.xlsx b/tests/gait/data/manual_gait_results_vcwt.xlsx
index 7ca3cb58..20e5ddbb 100644
Binary files a/tests/gait/data/manual_gait_results_vcwt.xlsx and b/tests/gait/data/manual_gait_results_vcwt.xlsx differ
diff --git a/tests/gait/data/manual_gait_turn_results_apcwt.xlsx b/tests/gait/data/manual_gait_turn_results_apcwt.xlsx
index e94faca4..24720e2c 100644
Binary files a/tests/gait/data/manual_gait_turn_results_apcwt.xlsx and b/tests/gait/data/manual_gait_turn_results_apcwt.xlsx differ
diff --git a/tests/gait/data/manual_gait_turn_results_vcwt.xlsx b/tests/gait/data/manual_gait_turn_results_vcwt.xlsx
index ad0e455b..f34fc6e5 100644
Binary files a/tests/gait/data/manual_gait_turn_results_vcwt.xlsx and b/tests/gait/data/manual_gait_turn_results_vcwt.xlsx differ
diff --git a/tests/gait/data/update_from_sheet_apcwt.py b/tests/gait/data/update_from_sheet_apcwt.py
index 2bd24e9d..320c1f04 100644
--- a/tests/gait/data/update_from_sheet_apcwt.py
+++ b/tests/gait/data/update_from_sheet_apcwt.py
@@ -56,8 +56,8 @@
     # get keys to update results for
     keys_turn = keys + ["Turn"]
 
-    tmp = {k: df.loc[:, k].values for k in keys}
-    tmp_turn = {k: df_turn.loc[:, k].values for k in keys_turn}
+    tmp = {k: df.loc[:, k].to_numpy(copy=True) for k in keys}
+    tmp_turn = {k: df_turn.loc[:, k].to_numpy(copy=True) for k in keys_turn}
 
     np.savez("gait_results_apcwt.npz", **tmp)
     np.savez("gait_results2_apcwt.npz", **tmp_turn)
diff --git a/tests/gait/data/update_from_sheet_vcwt.py b/tests/gait/data/update_from_sheet_vcwt.py
index 21ebb0b5..b077f4ca 100644
--- a/tests/gait/data/update_from_sheet_vcwt.py
+++ b/tests/gait/data/update_from_sheet_vcwt.py
@@ -88,12 +88,12 @@
     # old = np.load("gait_results_vcwt.npz")
     # old_turn = np.load("gait_results2_vcwt.npz")
 
-    tmp = {k: df.loc[:, k].values for k in keys}
+    tmp = {k: df.loc[:, k].to_numpy(copy=True) for k in keys}
     # tmp.update(
     #     {k: np.full(tmp["Bout N"].size, np.nan) for k in old.files if k not in tmp}
     # )
 
-    tmp_turn = {k: df_turn.loc[:, k].values for k in keys_turn}
+    tmp_turn = {k: df_turn.loc[:, k].to_numpy(copy=True) for k in keys_turn}
     # tmp_turn.update(
     #     {
     #         k: np.full(tmp_turn["Bout N"].size, np.nan)
diff --git a/tests/gait_old/data/gait_results.npz b/tests/gait_old/data/gait_results.npz
index 9b98496e..02fb1113 100644
Binary files a/tests/gait_old/data/gait_results.npz and b/tests/gait_old/data/gait_results.npz differ
diff --git a/tests/gait_old/data/gait_results2.npz b/tests/gait_old/data/gait_results2.npz
index 8cc53f06..8624e3ea 100644
Binary files a/tests/gait_old/data/gait_results2.npz and b/tests/gait_old/data/gait_results2.npz differ
diff --git a/tests/gait_old/data/manual_gait_results.xlsx b/tests/gait_old/data/manual_gait_results.xlsx
index d35df3bf..68768fad 100644
Binary files a/tests/gait_old/data/manual_gait_results.xlsx and b/tests/gait_old/data/manual_gait_results.xlsx differ
diff --git a/tests/gait_old/data/manual_gait_turn_results.xlsx b/tests/gait_old/data/manual_gait_turn_results.xlsx
index c39700a8..ca262388 100644
Binary files a/tests/gait_old/data/manual_gait_turn_results.xlsx and b/tests/gait_old/data/manual_gait_turn_results.xlsx differ
diff --git a/tests/gait_old/data/update_from_sheet.py b/tests/gait_old/data/update_from_sheet.py
index f1e45eb2..df3b2872 100644
--- a/tests/gait_old/data/update_from_sheet.py
+++ b/tests/gait_old/data/update_from_sheet.py
@@ -62,12 +62,12 @@
     old = np.load("gait_results.npz")
     old_turn = np.load("gait_results2.npz")
 
-    tmp = {k: df.loc[:, k].values for k in keys}
+    tmp = {k: df.loc[:, k].to_numpy(copy=True) for k in keys}
     # tmp.update(
     #     {k: np.full(tmp["Bout N"].size, np.nan) for k in old.files if k not in tmp}
     # )
 
-    tmp_turn = {k: df_turn.loc[:, k].values for k in keys_turn}
+    tmp_turn = {k: df_turn.loc[:, k].to_numpy(copy=True) for k in keys_turn}
     # tmp_turn.update(
     #     {
     #         k: np.full(tmp_turn["Bout N"].size, np.nan)
diff --git a/tests/io/test_read_base.py b/tests/io/test_read_base.py
index 16253391..c949209a 100644
--- a/tests/io/test_read_base.py
+++ b/tests/io/test_read_base.py
@@ -28,21 +28,19 @@ def test_in_pipeline(dummy_reader_class, dummy_process):
 def test_handle_naive_timestamps():
     # get a timestamp array of naive timestamps
     ts = (
-        date_range("2023-11-05 00:00:00", "2023-11-05 05:00:00", freq="0.1h").view(
-            "int64"
-        )
-        / 1e9
+        date_range(
+            "2023-11-05 00:00:00", "2023-11-05 05:00:00", freq="0.1h"
+        ).as_unit('s').view("int64")
     )
     # get the true timestamps - note that with the DST change there is an extra hour here
     # hence the earlier end time by 1 hour
     ts_true = (
         date_range(
-            "2023-11-05 00:00:00", "2023-11-05 04:00:00", freq="0.1h", tz="US/Eastern"
-        ).view("int64")
-        / 1e9
+            "2023-11-05 00:00:00", "2023-11-05 04:00:00", freq="0.1h", tz="America/New_York"
+        ).as_unit('s').view("int64")
     )
 
-    ts_pred = handle_naive_timestamps(ts, is_local=True, tz_name="US/Eastern")
+    ts_pred = handle_naive_timestamps(ts, is_local=True, tz_name="America/New_York")
 
     assert allclose(ts_pred - ts_true[0], ts_true - ts_true[0])
 
diff --git a/tests/io/test_read_csv.py b/tests/io/test_read_csv.py
index 705dbc6b..262a1b2f 100644
--- a/tests/io/test_read_csv.py
+++ b/tests/io/test_read_csv.py
@@ -24,10 +24,10 @@ def test_block_timestamps(self, dummy_csv_contents):
         )
 
         # get numpy arrays
-        time = raw["_datetime_"].astype(int).values / 1e9
+        time = raw["_datetime_"].dt.as_unit('s').astype(int).to_numpy(copy=True)
         data = {
-            "accel": raw[["ax", "ay", "az"]].values,
-            "temperature": raw["temperature"].values,
+            "accel": raw[["ax", "ay", "az"]].to_numpy(copy=True),
+            "temperature": raw["temperature"].to_numpy(copy=True),
         }
 
         comp_fs, time_rs, data_rs = rdr.handle_timestamp_inconsistency_np(
@@ -42,7 +42,7 @@ def test_block_timestamps(self, dummy_csv_contents):
             assert isnan(dstream).sum() == 0
 
         # trim a few samples off the last block and check we get a warning
-        time2 = raw["_datetime_"].astype(int).values[:-5] / 1e9
+        time2 = raw["_datetime_"].dt.as_unit('s').astype(int).values[:-5]
         data2 = {
             "accel": raw[["ax", "ay", "az"]].values[:-5],
             "temperature": raw["temperature"].values[:-5],
@@ -74,7 +74,7 @@ def test_unequal_blocks(self, dummy_csv_contents):
         raw.drop(index=range(13095, 14003), inplace=True)
         raw.reset_index(drop=True, inplace=True)
 
-        time = raw["_datetime_"].astype(int).values / 1e9
+        time = raw["_datetime_"].dt.as_unit('s').astype(int).values
         data = {
             "accel": raw[["ax", "ay", "az"]].values,
             "temperature": raw["temperature"].values,
diff --git a/tests/preprocessing/test_window.py b/tests/preprocessing/test_window.py
index e84e4203..ae46cb36 100644
--- a/tests/preprocessing/test_window.py
+++ b/tests/preprocessing/test_window.py
@@ -75,7 +75,7 @@ def test_dst(self):
             freq="s",
             tz="US/Eastern",
         )
-        time = dr.astype(int).values / 1e9
+        time = dr.as_unit('s').astype(int).values
 
         days = GetDayWindowIndices(bases=[0], periods=[24]).predict(
             time=time, fs=1.0, tz_name="US/Eastern"
@@ -93,7 +93,7 @@ def test_bad_guess(self):
             freq="s",
             tz="US/Eastern",
         )
-        time = dr.astype(int).values / 1e9
+        time = dr.as_unit('s').astype(int).values
 
         # intentionally making fs wrong here to force a bad guess at the end
         days = GetDayWindowIndices(bases=[0], periods=[24]).predict(