diff --git a/.github/workflows/build_publish.yml b/.github/workflows/build_publish.yml index adc3df8c..cbed893a 100644 --- a/.github/workflows/build_publish.yml +++ b/.github/workflows/build_publish.yml @@ -39,12 +39,13 @@ jobs: - name: Build wheels uses: pypa/cibuildwheel@v3.4.0 env: - CIBW_SKIP: pp* cp36-* cp37-* cp38-* cp39-* cp310-* *musllinux* + CIBW_SKIP: pp* cp38-* cp39-* cp310-* *musllinux* CIBW_ARCHS_MACOS: auto64 CIBW_ARCHS_LINUX: auto64 CIBW_ARCHS_WINDOWS: auto64 CIBW_BEFORE_BUILD_WINDOWS: "pip install delvewheel" CIBW_REPAIR_WHEEL_COMMAND_WINDOWS: "delvewheel repair -w {dest_dir} {wheel}" + MACOSX_DEPLOYMENT_TARGET: 15.0 with: output-dir: wheelhouse config-file: "{package}/pyproject.toml" @@ -69,6 +70,7 @@ jobs: path: dist/*.tar.gz upload_pypi: + if: ${{ github.event_name != 'workflow_dispatch' }} needs: [build_wheels, build_sdist] runs-on: ubuntu-latest steps: diff --git a/.readthedocs.yml b/.readthedocs.yml index ee75666c..44c816aa 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -25,9 +25,9 @@ python: extra_requirements: - docs build: - os: ubuntu-20.04 + os: ubuntu-lts-latest tools: - python: "3.10" + python: "3.13" jobs: pre_install: - echo "\n\n[tool.meson-python.args]\nsetup = [\"-Dis_not_rtd_build=false\"]" >> pyproject.toml diff --git a/meson.build b/meson.build index 120f67bc..51043d2d 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project( 'scikit-digital-health', 'c', - version: '0.17.11', + version: '0.17.12', license: 'MIT', meson_version: '>=1.1', ) diff --git a/src/skdh/activity/core.py b/src/skdh/activity/core.py index a668a145..c6ced038 100644 --- a/src/skdh/activity/core.py +++ b/src/skdh/activity/core.py @@ -209,6 +209,10 @@ def __init__( self.wake_endpoints = [ ept.IntensityGradient(state="wake"), ept.MaxAcceleration(self.max_acc_lens, state="wake"), + # M5hr, M10hr + ept.MaxAcceleration([300, 600], required_points=0.5, state="wake"), + # L5hr + ept.MinAcceleration([300], required_points=0.5, state="wake"), ] self.wake_endpoints += [ ept.TotalIntensityTime(lvl, self.wlen, self.cutpoints, state="wake") diff --git a/src/skdh/activity/endpoints.py b/src/skdh/activity/endpoints.py index bb046eb7..50f6d80d 100644 --- a/src/skdh/activity/endpoints.py +++ b/src/skdh/activity/endpoints.py @@ -11,7 +11,9 @@ zeros, full, max, + min, nanmax, + nanmin, histogram, log, nan, @@ -43,6 +45,7 @@ "ActivityEndpoint", "IntensityGradient", "MaxAcceleration", + "MinAcceleration", "TotalIntensityTime", "BoutIntensityTime", "FragmentationEndpoints", @@ -394,17 +397,23 @@ class MaxAcceleration(ActivityEndpoint): ---------- window_lengths : {list, int} List of window lengths, or a single window length. + required_points : int, float, optional + Number of points required in the window to compute the endpoint. If float, + should be between 0 and 1 and will be interpreted as a fraction of the total + number of points in the window. Default is 1.0, meaning all points in the + window must be present. state : {'wake', 'sleep} State during which the endpoint is being computed. """ - def __init__(self, window_lengths, state="wake"): + def __init__(self, window_lengths, required_points=1.0, state="wake"): if isinstance(window_lengths, int): window_lengths = [window_lengths] super().__init__([f"max acc {i}min [g]" for i in window_lengths], state) self.wlens = window_lengths + self.required_points = required_points def predict( self, @@ -445,7 +454,7 @@ def predict( # skipping more samples would introduce bias by random chance of # where the windows start and stop try: - tmp_max = max(moving_mean(accel_metric, n, 1)) + tmp_max = max(moving_mean(accel_metric, n, 1, req_points=self.required_points)) except ValueError: return # if the window length is too long for this block of data @@ -453,6 +462,79 @@ def predict( results[name][i] = nanmax([tmp_max, results[name][i]]) +class MinAcceleration(ActivityEndpoint): + """ + Compute the minimum acceleration over windows of the specified length. + + Parameters + ---------- + window_lengths : {list, int} + List of window lengths, or a single window length. + required_points : int, float, optional + Number of points required in the window to compute the endpoint. If float, + should be between 0 and 1 and will be interpreted as a fraction of the total + number of points in the window. Default is 1.0, meaning all points in the + window must be present. + state : {'wake', 'sleep} + State during which the endpoint is being computed. + """ + + def __init__(self, window_lengths, required_points=1.0, state="wake"): + if isinstance(window_lengths, int): + window_lengths = [window_lengths] + + super().__init__([f"min acc {i}min [g]" for i in window_lengths], state) + + self.wlens = window_lengths + self.required_points = required_points + + def predict( + self, + results, + i, + accel_metric, + accel_metric_60, + epoch_s, + epochs_per_min, + **kwargs, + ): + """ + Compute the minimum acceleration during this set of data, and compare it + to the previous smallest detected value. + + Parameters + ---------- + results : dict + Dictionary containing the initialized results arrays. Keys in `results` + are taken from the names of endpoints. + i : int + Index of the day, used to index into individual result arrays, e.g. + `results[self.name][i] = 5.0` + accel_metric : numpy.ndarray + Computed acceleration metric (e.g. ENMO). + accel_metric_60 : numpy.ndarray + Computed acceleration metric for a 60 second window. + epoch_s : int + Duration in seconds of each sample of `accel_metric`. + epochs_per_min : int + Number of epochs per minute. + """ + super(MinAcceleration, self).predict() + + for wlen, name in zip(self.wlens, self.name): + n = wlen * epochs_per_min + # skip 1 sample because we want the window with the smallest acceleration + # skipping more samples would introduce bias by random chance of + # where the windows start and stop + try: + tmp_min = min(moving_mean(accel_metric, n, 1, req_points=self.required_points)) + except ValueError: + return # if the window length is too long for this block of data + + # check that we don't have a smaller result already for this day + results[name][i] = nanmin([tmp_min, results[name][i]]) + + class TotalIntensityTime(ActivityEndpoint): """ Compute the total time spent in an intensity level. diff --git a/src/skdh/features/lib/extensions/_utility.c b/src/skdh/features/lib/extensions/_utility.c index bcca7d2a..9d8dc062 100644 --- a/src/skdh/features/lib/extensions/_utility.c +++ b/src/skdh/features/lib/extensions/_utility.c @@ -7,18 +7,18 @@ #include #include -extern void mean_sd_1d(long *, double *, double *, double *); -extern void unique(long *, double *, double *, long *, long *); -extern void gmean(long *, double *, double *); -extern void embed_sort(long *, long *, double *, long *, long *, long *); -extern void hist(long *, double *, long *, double *, double *, long *); -extern void histogram(long *, long *, double *, double *, long *); +extern void mean_sd_1d(Py_ssize_t *, double *, double *, double *); +extern void unique(Py_ssize_t *, double *, double *, long *, long *); +extern void gmean(Py_ssize_t *, double *, double *); +extern void embed_sort(Py_ssize_t *, long *, double *, long *, long *, long *); +extern void hist(Py_ssize_t *, double *, Py_ssize_t *, double *, double *, long *); +extern void histogram(Py_ssize_t *, Py_ssize_t *, double *, double *, long *); extern void insertion_sort_2d(long *, long *, double *, long *); extern void quick_argsort_(long *, double *, long *); extern void quick_sort_(long *, double *); -extern void f_rfft(long *, double *, long *, double *); +extern void f_rfft(Py_ssize_t *, double *, long *, double *); extern void destroy_plan(void); @@ -194,10 +194,10 @@ PyObject * cf_embed_sort(PyObject *NPY_UNUSED(self), PyObject *args){ PyObject * cf_hist(PyObject *NPY_UNUSED(self), PyObject *args){ PyObject *x_; - long ncells = 1; + npy_intp ncells = 1; double min_val = 0., max_val = 100.; - if (!PyArg_ParseTuple(args, "Oldd:cf_hist", &x_, &ncells, &min_val, &max_val)) return NULL; + if (!PyArg_ParseTuple(args, "Ondd:cf_hist", &x_, &ncells, &min_val, &max_val)) return NULL; PyArrayObject *data = (PyArrayObject *)PyArray_FromAny( x_, PyArray_DescrFromType(NPY_DOUBLE), 1, 0, @@ -261,8 +261,8 @@ PyObject * cf_histogram(PyObject *NPY_UNUSED(self), PyObject *args){ npy_intp n_elem = PyArray_SIZE(data); - long n3 = 3; - long ncells = (long)ceil(sqrt(n_elem)); + npy_intp n3 = 3; + npy_intp ncells = (npy_intp)ceil(sqrt(n_elem)); PyArrayObject *d = (PyArrayObject *)PyArray_EMPTY(1, &n3, NPY_DOUBLE, 0); PyArrayObject *counts = (PyArrayObject *)PyArray_ZEROS(1, &ncells, NPY_LONG, 0); diff --git a/src/skdh/features/lib/extensions/entropy.c b/src/skdh/features/lib/extensions/entropy.c index 51fa9715..26134ef2 100644 --- a/src/skdh/features/lib/extensions/entropy.c +++ b/src/skdh/features/lib/extensions/entropy.c @@ -6,9 +6,9 @@ #include #include -extern void signal_entropy_1d(long *, double *, double *); -extern void sample_entropy_1d(long *, double *, long *, double *, double *); -extern void permutation_entropy_1d(long *, double *, long *, long *, int *, double *); +extern void signal_entropy_1d(Py_ssize_t *, double *, double *); +extern void sample_entropy_1d(Py_ssize_t *, double *, long *, double *, double *); +extern void permutation_entropy_1d(Py_ssize_t *, double *, long *, long *, int *, double *); PyObject * signal_entropy(PyObject *NPY_UNUSED(self), PyObject *args){ @@ -48,7 +48,7 @@ PyObject * signal_entropy(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -107,7 +107,7 @@ PyObject * sample_entropy(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -168,7 +168,7 @@ PyObject * permutation_entropy(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ diff --git a/src/skdh/features/lib/extensions/f_rfft.f95 b/src/skdh/features/lib/extensions/f_rfft.f95 index 42d4fd27..046a9257 100644 --- a/src/skdh/features/lib/extensions/f_rfft.f95 +++ b/src/skdh/features/lib/extensions/f_rfft.f95 @@ -6,7 +6,8 @@ subroutine f_rfft(n, x, nfft, F) bind(C, name="f_rfft") use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long), intent(in) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: F(2 * nfft + 2) ! local diff --git a/src/skdh/features/lib/extensions/ffeatures.f95 b/src/skdh/features/lib/extensions/ffeatures.f95 index 4729408e..d598c7af 100644 --- a/src/skdh/features/lib/extensions/ffeatures.f95 +++ b/src/skdh/features/lib/extensions/ffeatures.f95 @@ -7,7 +7,7 @@ ! Compute the autocorrelation of a signal with the specified lag ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute autocorrelation for ! lag : integer(long), lag for the autocorrelation, in samples ! normalize : integer(int), normalize the autocorrelation @@ -19,30 +19,33 @@ subroutine autocorr_1d(n, x, lag, normalize, ac) bind(C, name="autocorr_1d") use, intrinsic :: iso_c_binding use utility, only : mean_sd_1d implicit none - integer(c_long), intent(in) :: n, lag + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: lag real(c_double), intent(in) :: x(n) integer(c_int), intent(in) :: normalize real(c_double), intent(out) :: ac ! local - integer(c_long) :: i + integer(c_size_t) :: i, st_lag real(c_double) :: mean1, mean2, std1, std2 + + st_lag = int(lag, c_size_t) ac = 0._c_double if (normalize == 1_c_int) then - call mean_sd_1d(n - lag, x(:n - lag), mean1, std1) - call mean_sd_1d(n - lag, x(lag + 1:), mean2, std2) + call mean_sd_1d(n - st_lag, x(:n - st_lag), mean1, std1) + call mean_sd_1d(n - st_lag, x(st_lag + 1:), mean2, std2) - do i=1_c_long, n - lag - ac = ac + (x(i) - mean1) * (x(i + lag) - mean2) + do i=1_c_size_t, n - st_lag + ac = ac + (x(i) - mean1) * (x(i + st_lag) - mean2) end do ac = ac / ((n - lag - 1) * std1 * std2) else - call mean_sd_1d(n - lag, x(:n - lag), mean1, std1) - call mean_sd_1d(n - lag, x(lag + 1:), mean2, std2) + call mean_sd_1d(n - st_lag, x(:n - st_lag), mean1, std1) + call mean_sd_1d(n - st_lag, x(st_lag + 1:), mean2, std2) - do i=1_c_long, n - lag - ac = ac + x(i) * x(i + lag) + do i=1_c_size_t, n - st_lag + ac = ac + x(i) * x(i + st_lag) end do ac = ac / (std1 * std2) end if @@ -54,7 +57,7 @@ subroutine autocorr_1d(n, x, lag, normalize, ac) bind(C, name="autocorr_1d") ! Compute the complexity invariant distance metric for a signal ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute over ! normalize : integer(int), normalize the complexity invariant distance ! @@ -65,12 +68,12 @@ subroutine cid_1d(n, x, normalize, cid) bind(C, name='cid_1d') use, intrinsic :: iso_c_binding use utility, only : mean_sd_1d implicit none - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n) integer(c_int), intent(in) :: normalize real(c_double), intent(out) :: cid ! local - integer(c_long) :: i + integer(c_size_t) :: i real(c_double) :: mean, sd cid = 0._c_double; @@ -104,7 +107,8 @@ subroutine permutation_entropy_1d(n, x, order, delay, normalize, pe) bind(C, nam use, intrinsic :: iso_c_binding use utility, only : embed_sort, unique implicit none - integer(c_long), intent(in) :: n, order, delay + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: order, delay real(c_double), intent(in) :: x(n) integer(c_int), intent(in) :: normalize real(c_double), intent(out) :: pe @@ -116,7 +120,7 @@ subroutine permutation_entropy_1d(n, x, order, delay, normalize, pe) bind(C, nam real(c_double) :: unq_counts(n - (order - 1) * delay) real(c_double), parameter :: log2 = dlog(2._c_double) - nsi = n - (order - 1) * delay + nsi = int(n, c_long) - (order - 1) * delay call embed_sort(n, nsi, x, order, delay, sorted_idx) @@ -126,7 +130,7 @@ subroutine permutation_entropy_1d(n, x, order, delay, normalize, pe) bind(C, nam end do unq_counts = 0._c_double - call unique(nsi, hashval, unq_vals, unq_counts, n_unq) + call unique(int(nsi, c_size_t), hashval, unq_vals, unq_counts, n_unq) unq_counts = unq_counts / sum(unq_counts) @@ -143,7 +147,7 @@ subroutine permutation_entropy_1d(n, x, order, delay, normalize, pe) bind(C, nam ! Compute the sample entropy of a signal ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute sample entropy on ! L : integer(long), length of sets to compare ! r : real(double), maximum set distance @@ -155,7 +159,8 @@ subroutine sample_entropy_1d(n, x, L, r, samp_ent) bind(C, name="sample_entropy_ use, intrinsic :: iso_c_binding use, intrinsic :: iso_fortran_env, only: stdout=>output_unit implicit none - integer(c_long), intent(in) :: n, L + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: L real(c_double), intent(in) :: x(n), r real(c_double), intent(out) :: samp_ent ! local @@ -199,7 +204,7 @@ subroutine sample_entropy_1d(n, x, L, r, samp_ent) bind(C, name="sample_entropy_ ! Compute the signal entropy of a 1d signal ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute signal entropy for ! ! Output @@ -209,13 +214,13 @@ subroutine signal_entropy_1d(n, x, ent) bind(C, name="signal_entropy_1d") use, intrinsic :: iso_c_binding use utility, only : histogram, mean_sd_1d implicit none - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(out) :: x(n) real(c_double), intent(out) :: ent ! local real(c_double) :: d(3), nbias, cnt real(c_double) :: estimate, stdev, mean, logf - integer(c_long) :: i, sqrt_n + integer(c_size_t) :: i, sqrt_n integer(c_long) :: h(ceiling(sqrt(real(n)))) sqrt_n = ceiling(sqrt(real(n, c_double))) @@ -237,7 +242,7 @@ subroutine signal_entropy_1d(n, x, ent) bind(C, name="signal_entropy_1d") cnt = 0._c_double estimate = 0._c_double - do i = 1, int(d(3), c_long) + do i = 1, int(d(3), c_size_t) if (h(i) > 0) then logf = log(real(h(i), c_double)) else @@ -260,7 +265,7 @@ subroutine signal_entropy_1d(n, x, ent) bind(C, name="signal_entropy_1d") ! Compute the dominant frequency of a signal, in the specified range ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute signal entropy for ! nfft : integer(long), number of points to use in the FFT computation ! fs : real(double), sampling frequency in Hz @@ -274,7 +279,8 @@ subroutine dominant_freq_1d(n, x, fs, nfft, low_cut, hi_cut, df) bind(C, name="d use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n), low_cut, hi_cut, fs real(c_double), intent(out) :: df ! local @@ -312,7 +318,7 @@ subroutine dominant_freq_1d(n, x, fs, nfft, low_cut, hi_cut, df) bind(C, name="d ! Compute the dominant frequency value (spectral power) of a signal, in the specified range ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute signal entropy for ! nfft : integer(long), number of points to use in the FFT computation ! fs : real(double), sampling frequency in Hz @@ -326,7 +332,8 @@ subroutine dominant_freq_value_1d(n, x, fs, nfft, low_cut, hi_cut, dfval) bind(C use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n), low_cut, hi_cut, fs real(c_double), intent(out) :: dfval ! local @@ -361,7 +368,7 @@ subroutine dominant_freq_value_1d(n, x, fs, nfft, low_cut, hi_cut, dfval) bind(C ! dominant frequency ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute signal entropy for ! nfft : integer(long), number of points to use in the FFT computation ! fs : real(double), sampling frequency in Hz @@ -375,7 +382,8 @@ subroutine power_spectral_sum_1d(n, x, fs, nfft, low_cut, hi_cut, pss) bind(C, n use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long), intent(in) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n), fs, low_cut, hi_cut real(c_double), intent(out) :: pss ! local @@ -421,7 +429,7 @@ subroutine power_spectral_sum_1d(n, x, fs, nfft, low_cut, hi_cut, pss) bind(C, n ! Compute the sum of the spectral power in the range specified. ! ! Input -! n : integer(long) +! n : integer(c_size_t) ! x(n) : real(double), array to compute signal entropy for ! nfft : integer(long), number of points to use in the FFT computation ! fs : real(double), sampling frequency in Hz @@ -438,13 +446,14 @@ subroutine range_power_sum_1d(n, x, fs, nfft, low_cut, hi_cut, demean, use_mod, use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long), intent(in) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n), fs, low_cut, hi_cut integer(c_int), intent(in) :: demean, use_mod, normalize real(c_double), intent(out) :: rps ! local real(c_double), parameter :: log2 = log(2._c_double) - integer(c_long) :: i, ihcut, ilcut, ier, istart + integer(c_long) :: ihcut, ilcut, ier, istart real(c_double) :: sp_norm(nfft + 1) real(c_double) :: sp_hat(2 * nfft + 2), y(2 * nfft) @@ -489,7 +498,7 @@ subroutine range_power_sum_1d(n, x, fs, nfft, low_cut, hi_cut, demean, use_mod, ! Compute the spectral entropy of the specified frequency range ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute signal entropy for ! nfft : integer(long), number of points to use in the FFT computation ! fs : real(double), sampling frequency in Hz @@ -503,7 +512,8 @@ subroutine spectral_entropy_1d(n, x, fs, nfft, low_cut, hi_cut, sEnt) bind(C, na use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long), intent(in) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n), low_cut, hi_cut, fs real(c_double), intent(out) :: sEnt ! local @@ -542,7 +552,7 @@ subroutine spectral_entropy_1d(n, x, fs, nfft, low_cut, hi_cut, sEnt) bind(C, na ! Compute the spectral flatness of the specified frequency range ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute signal entropy for ! nfft : integer(long), number of points to use in the FFT computation ! fs : real(double), sampling frequency in Hz @@ -557,7 +567,8 @@ subroutine spectral_flatness_1d(n, x, fs, nfft, low_cut, hi_cut, sFlat) bind(C, use real_fft, only : execute_real_forward use utility, only : gmean implicit none - integer(c_long), intent(in) :: n, nfft + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: nfft real(c_double), intent(in) :: x(n), low_cut, hi_cut, fs real(c_double), intent(out) :: sFlat ! local @@ -583,7 +594,7 @@ subroutine spectral_flatness_1d(n, x, fs, nfft, low_cut, hi_cut, sFlat) bind(C, sp_norm = sp_norm / sum(sp_norm(ilcut:ihcut)) + 1.d-10 mean = sum(sp_norm(ilcut:ihcut)) / (ihcut - ilcut + 1) - call gmean(ihcut - ilcut + 1, sp_norm(ilcut:ihcut), sFlat) + call gmean(int(ihcut - ilcut + 1, c_size_t), sp_norm(ilcut:ihcut), sFlat) sFlat = 10._c_double * log(sFlat / mean) / log(10._c_double) end subroutine @@ -593,7 +604,7 @@ subroutine spectral_flatness_1d(n, x, fs, nfft, low_cut, hi_cut, sFlat) bind(C, ! Compute the jerk metric for a 1d signal ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute for ! fs : real(double), sampling frequency, in Hz ! @@ -603,11 +614,11 @@ subroutine spectral_flatness_1d(n, x, fs, nfft, low_cut, hi_cut, sFlat) bind(C, subroutine jerk_1d(n, x, fs, jerk) bind(C, name="jerk_1d") use, intrinsic :: iso_c_binding implicit none - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n), fs real(c_double), intent(out) :: jerk ! local - integer(c_long) :: i + integer(c_size_t) :: i real(c_double) :: amp, jsum, xold jsum = 0._c_double @@ -629,7 +640,7 @@ subroutine jerk_1d(n, x, fs, jerk) bind(C, name="jerk_1d") ! Compute the jerk metric, but dimensionless, for a signal ! ! Input -! n : integer(long), axis dimension +! n : integer(size_t), axis dimension ! x(n) : real(double), array to compute for ! stype : integer(long), type of signal input. 1=velocity, 2=acceleration, 3=jerk ! @@ -639,11 +650,12 @@ subroutine jerk_1d(n, x, fs, jerk) bind(C, name="jerk_1d") subroutine dimensionless_jerk_1d(n, x, stype, djerk) bind(C, name="dimensionless_jerk_1d") use, intrinsic :: iso_c_binding implicit none - integer(c_long), intent(in) :: n, stype + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: stype real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: djerk ! local - integer(c_long) :: i + integer(c_size_t) :: i real(c_double) :: amp, jsum jsum = 0._c_double @@ -685,7 +697,7 @@ subroutine dimensionless_jerk_1d(n, x, stype, djerk) bind(C, name="dimensionless ! Compute the percentage of samples that lie within the specified range ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute for ! xmin : real(double), range min value ! xmax : real(double), range max value @@ -696,11 +708,11 @@ subroutine dimensionless_jerk_1d(n, x, stype, djerk) bind(C, name="dimensionless subroutine range_count_1d(n, x, xmin, xmax, rcp) bind(C, name="range_count_1d") use, intrinsic :: iso_c_binding implicit none - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n), xmin, xmax real(c_double), intent(out) :: rcp ! local - integer(c_long) :: i + integer(c_size_t) :: i rcp = 0._c_double do i=1, n @@ -718,7 +730,7 @@ subroutine range_count_1d(n, x, xmin, xmax, rcp) bind(C, name="range_count_1d") ! Compute the percentage of samples farther away than r * SD from the mean ! ! Input -! n : integer(long) +! n : integer(size_t) ! x(n) : real(double), array to compute for ! r : real(double), factor to multiply SD by ! @@ -729,7 +741,7 @@ subroutine ratio_beyond_r_sigma_1d(n, x, r, rbrs) bind(C, name="ratio_beyond_r_s use, intrinsic :: iso_c_binding use utility, only : mean_sd_1d implicit none - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n), r real(c_double), intent(out) :: rbrs ! local @@ -748,7 +760,7 @@ subroutine ratio_beyond_r_sigma_1d(n, x, r, rbrs) bind(C, name="ratio_beyond_r_s ! Compute the linear regression and the slope ! ! Input -! n : integer(long) +! n : integer(size_t) ! y(n) : real(double), array to compute for ! fs : real(double), sampling frequency for the array ! @@ -758,11 +770,11 @@ subroutine ratio_beyond_r_sigma_1d(n, x, r, rbrs) bind(C, name="ratio_beyond_r_s subroutine linear_regression_1d(n, y, fs, slope) bind(C, name="linear_regression_1d") use, intrinsic :: iso_c_binding implicit none - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: y(n), fs real(c_double), intent(out) :: slope ! local - integer(c_long) :: i + integer(c_size_t) :: i real(c_double) :: ssxm, ssxym, ky, Ex, Ey, Exy ! The variance of a time series with N elements and sampling frequency fs @@ -788,7 +800,7 @@ subroutine linear_regression_1d(n, y, fs, slope) bind(C, name="linear_regression Ey = 0._c_double Exy = 0._c_double - do i=2, n + do i=2_c_size_t, n Ey = Ey + (y(i) - ky) Exy = Exy + ((i-1) / fs) * (y(i) - ky) end do @@ -803,7 +815,7 @@ subroutine linear_regression_1d(n, y, fs, slope) bind(C, name="linear_regression ! Compute the spectral arc length measure of smoothness ! ! Input -! n : integer(long), axis dimension +! n : integer(size_t), axis dimension ! x(n) : real(double), array to compute for ! fs : real(double), sampling frequency in Hz ! padlevel : integer(long), amount of zero-padding for the FFT @@ -818,11 +830,13 @@ subroutine sparc_1d(n, x, fs, padlevel, fc, amp_thresh, sal) bind(C, name="sparc use, intrinsic :: iso_c_binding use real_fft, only : execute_real_forward implicit none - integer(c_long), intent(in) :: n, padlevel + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: padlevel real(c_double), intent(in) :: x(n), fs, fc, amp_thresh real(c_double), intent(out) :: sal ! local - integer(c_long) :: nfft, i, ixf, inxi, inxf, ier + integer(c_size_t) :: nfft, i + integer(c_long) :: ixf, inxi, inxf, ier real(c_double) :: freq_factor real(c_double) :: Mf(2**(ceiling(log(real(n))/log(2.0)) + padlevel-1)+1) real(c_double) :: y(2**(ceiling(log(real(n))/log(2.0)) + padlevel)) @@ -842,7 +856,7 @@ subroutine sparc_1d(n, x, fs, padlevel, fc, amp_thresh, sal) bind(C, name="sparc sp_hat = 0._c_double y = 0._c_double y(:n) = x - call execute_real_forward(nfft, y, 1._c_double, sp_hat, ier) + call execute_real_forward(int(nfft, c_long), y, 1._c_double, sp_hat, ier) if (ier /= 0_c_long) return ! normalize the FFT response diff --git a/src/skdh/features/lib/extensions/frequency.c b/src/skdh/features/lib/extensions/frequency.c index f62d01cc..3c73e2c6 100644 --- a/src/skdh/features/lib/extensions/frequency.c +++ b/src/skdh/features/lib/extensions/frequency.c @@ -7,12 +7,12 @@ #include #include -extern void dominant_freq_1d(long *, double *, double *, long *, double *, double *, double *); -extern void dominant_freq_value_1d(long *, double *, double *, long *, double *, double *, double *); -extern void power_spectral_sum_1d(long *, double *, double *, long *, double *, double *, double *); -extern void range_power_sum_1d(long *, double *, double *, long *, double *, double *, int *, int *, int *, double *); -extern void spectral_entropy_1d(long *, double *, double *, long *, double *, double *, double *); -extern void spectral_flatness_1d(long *, double *, double *, long *, double *, double *, double *); +extern void dominant_freq_1d(Py_ssize_t *, double *, double *, long *, double *, double *, double *); +extern void dominant_freq_value_1d(Py_ssize_t *, double *, double *, long *, double *, double *, double *); +extern void power_spectral_sum_1d(Py_ssize_t *, double *, double *, long *, double *, double *, double *); +extern void range_power_sum_1d(Py_ssize_t *, double *, double *, long *, double *, double *, int *, int *, int *, double *); +extern void spectral_entropy_1d(Py_ssize_t *, double *, double *, long *, double *, double *, double *); +extern void spectral_flatness_1d(Py_ssize_t *, double *, double *, long *, double *, double *, double *); extern void destroy_plan(void); @@ -66,7 +66,7 @@ PyObject * dominant_frequency(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -140,7 +140,7 @@ PyObject * dominant_frequency_value(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -215,7 +215,7 @@ PyObject * power_spectral_sum(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -298,7 +298,7 @@ PyObject * range_power_sum(PyObject *NPY_UNUSED(self), PyObject *args) double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i) @@ -376,7 +376,7 @@ PyObject * spectral_entropy(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -451,7 +451,7 @@ PyObject * spectral_flatness(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ diff --git a/src/skdh/features/lib/extensions/misc_features.c b/src/skdh/features/lib/extensions/misc_features.c index 81e00a18..9c80be45 100644 --- a/src/skdh/features/lib/extensions/misc_features.c +++ b/src/skdh/features/lib/extensions/misc_features.c @@ -6,9 +6,9 @@ #include #include -extern void cid_1d(long *, double *, int *, double *); -extern void range_count_1d(long *, double *, double *, double *, double *); -extern void ratio_beyond_r_sigma_1d(long *, double *, double *, double *); +extern void cid_1d(Py_ssize_t *, double *, int *, double *); +extern void range_count_1d(Py_ssize_t *, double *, double *, double *, double *); +extern void ratio_beyond_r_sigma_1d(Py_ssize_t *, double *, double *, double *); PyObject * complexity_invariant_distance(PyObject *NPY_UNUSED(self), PyObject *args){ PyObject *x_; @@ -53,7 +53,7 @@ PyObject * complexity_invariant_distance(PyObject *NPY_UNUSED(self), PyObject *a double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim - 1]; + npy_intp stride = ddims[ndim - 1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -111,7 +111,7 @@ PyObject * range_count(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -169,7 +169,7 @@ PyObject * ratio_beyond_r_sigma(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ diff --git a/src/skdh/features/lib/extensions/real_fft.f95 b/src/skdh/features/lib/extensions/real_fft.f95 index 25e83db9..061b3be1 100644 --- a/src/skdh/features/lib/extensions/real_fft.f95 +++ b/src/skdh/features/lib/extensions/real_fft.f95 @@ -676,7 +676,7 @@ subroutine rfftp_factorize(ier) plan%fct(nfct)%fct = tmp end if - maxl = int(sqrt(real( length)), 8) + 1 + maxl = int(sqrt(real( length)), 8) + 1_c_long divisor = 3_c_long do while (( length > 1) .AND. (divisor < maxl)) @@ -692,7 +692,7 @@ subroutine rfftp_factorize(ier) length = length / divisor end do - maxl = int(sqrt(real( length)), 8) + 1 + maxl = int(sqrt(real( length)), 8) + 1_c_long end if divisor = divisor + 2 diff --git a/src/skdh/features/lib/extensions/smoothness.c b/src/skdh/features/lib/extensions/smoothness.c index db6a42fe..339b8f71 100644 --- a/src/skdh/features/lib/extensions/smoothness.c +++ b/src/skdh/features/lib/extensions/smoothness.c @@ -7,9 +7,9 @@ #include -extern void jerk_1d(long *, double *, double *, double *); -extern void dimensionless_jerk_1d(long *, double *, long *, double *); -extern void sparc_1d(long *, double *, double *, long *, double *, double *, double *); +extern void jerk_1d(Py_ssize_t *, double *, double *, double *); +extern void dimensionless_jerk_1d(Py_ssize_t *, double *, long *, double *); +extern void sparc_1d(Py_ssize_t *, double *, double *, long *, double *, double *, double *); extern void destroy_plan(void); PyObject * jerk_metric(PyObject *NPY_UNUSED(self), PyObject *args){ @@ -50,7 +50,7 @@ PyObject * jerk_metric(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -108,7 +108,7 @@ PyObject * dimensionless_jerk_metric(PyObject *NPY_UNUSED(self), PyObject *args) double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -167,7 +167,7 @@ PyObject * SPARC(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ diff --git a/src/skdh/features/lib/extensions/sort.f95 b/src/skdh/features/lib/extensions/sort.f95 index 5ab0abd4..b3751342 100644 --- a/src/skdh/features/lib/extensions/sort.f95 +++ b/src/skdh/features/lib/extensions/sort.f95 @@ -74,11 +74,11 @@ subroutine insertion_sort_2d(m, order, x, idx) bind(C, name="insertion_sort_2d") ! idx(n) : integer(8), index of values to be rearranged (inplace) to match the sorting of x ! -------------------------------------------------------------------- recursive subroutine quick_argsort_(n, x, idx) bind(C, name="quick_argsort_") - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(inout) :: x(n) integer(c_long), intent(inout) :: idx(n) ! local - integer(c_long) :: i, j, k, l, g, itemp, ip1, ip2 + integer(c_size_t) :: i, j, k, l, g, itemp, ip1, ip2 real(c_double) :: temp, p1, p2 ! use insertion sort on small arrays @@ -202,10 +202,10 @@ recursive subroutine quick_argsort_(n, x, idx) bind(C, name="quick_argsort_") ! x(n) : real(8), array of values to sort (inplace) ! -------------------------------------------------------------------- recursive subroutine quick_sort_(n, x) bind(C, name="quick_sort_") - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(inout) :: x(n) !f2py intent(hide) :: n - integer(c_long) :: i, j, k, l, g + integer(c_size_t) :: i, j, k, l, g real(c_double) :: temp, p1, p2 ! use insertion sort on small arrays diff --git a/src/skdh/features/lib/extensions/statistics.c b/src/skdh/features/lib/extensions/statistics.c index 8247457e..cdce2741 100644 --- a/src/skdh/features/lib/extensions/statistics.c +++ b/src/skdh/features/lib/extensions/statistics.c @@ -7,8 +7,8 @@ #include -extern void autocorr_1d(long *, double *, long *, int *, double *); -extern void linear_regression_1d(long *, double *, double *, double *); +extern void autocorr_1d(Py_ssize_t *, double *, long *, int *, double *); +extern void linear_regression_1d(Py_ssize_t *, double *, double *, double *); PyObject * autocorrelation(PyObject *NPY_UNUSED(self), PyObject *args){ @@ -55,7 +55,7 @@ PyObject * autocorrelation(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ @@ -113,7 +113,7 @@ PyObject * linear_regression(PyObject *NPY_UNUSED(self), PyObject *args){ double *dptr = (double *)PyArray_DATA(data); double *rptr = (double *)PyArray_DATA(res); - long stride = ddims[ndim-1]; + npy_intp stride = ddims[ndim-1]; int nrepeats = PyArray_SIZE(data) / stride; for (int i = 0; i < nrepeats; ++i){ diff --git a/src/skdh/features/lib/extensions/utility.f95 b/src/skdh/features/lib/extensions/utility.f95 index 06ae762b..d82b6d7e 100644 --- a/src/skdh/features/lib/extensions/utility.f95 +++ b/src/skdh/features/lib/extensions/utility.f95 @@ -18,12 +18,12 @@ module utility ! sd : real(double), sample standard deviation of x ! -------------------------------------------------------------------- subroutine mean_sd_1d(n, x, mn, sd) bind(C) - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: mn, sd ! local real(c_double) :: Ex, Ex2 - integer(c_long) :: i + integer(c_size_t) :: i Ex2 = 0._c_double mn = 0._c_double @@ -55,12 +55,12 @@ subroutine mean_sd_1d(n, x, mn, sd) bind(C) ! j : integer(long), number of unique values in x ! -------------------------------------------------------------------- subroutine unique(n, x, uniq, counts, j) bind(C) - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: uniq(n), counts(n) integer(c_long), intent(out) :: j !f2py intent(hide) :: n - integer(c_long) :: i + integer(c_size_t) :: i real(c_double) :: y(n) y = x @@ -72,7 +72,7 @@ subroutine unique(n, x, uniq, counts, j) bind(C) counts(1) = 1._c_double j = 2_c_long - do i=2_c_long, n + do i=2_c_size_t, n if (y(i) .NE. y(i-1)) then uniq(j) = y(i) counts(j) = counts(j) + 1 @@ -90,25 +90,25 @@ subroutine unique(n, x, uniq, counts, j) bind(C) ! Compute the geometric mean of a series ! ! In - ! n : integer(long), number of samples in x + ! n : integer(size_t), number of samples in x ! x(n) : real(double), array of values of length n ! ! Out ! gm : real(double), geometric mean of the series x ! -------------------------------------------------------------------- subroutine gmean(n, x, gm) bind(C) - integer(c_long), intent(in) :: n + integer(c_size_t), intent(in) :: n real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: gm !f2py intent(hide) :: n real(c_double) :: logsum, prod real(c_double), parameter :: large=1.d64, small=1.d-64 - integer(c_long) :: i + integer(c_size_t) :: i logsum = 0._c_double prod = 1._c_double - do i=1, n + do i=1_c_size_t, n prod = prod * x(i) if ((prod > large) .OR. (prod < small)) then logsum = logsum + log(prod) @@ -132,7 +132,8 @@ subroutine gmean(n, x, gm) bind(C) ! res(m, order) : integer(long), sorted embedding vectors ! -------------------------------------------------------------------- subroutine embed_sort(n, m, x, order, delay, res) bind(C, name="embed_sort") - integer(c_long), intent(in) :: n, m, order, delay + integer(c_size_t), intent(in) :: n + integer(c_long), intent(in) :: m, order, delay real(c_double), intent(in) :: x(n) integer(c_long), intent(out) :: res(m, order) ! local @@ -165,11 +166,11 @@ subroutine embed_sort(n, m, x, order, delay, res) bind(C, name="embed_sort") ! counts : integer(8), counts in each bin ! -------------------------------------------------------------------- subroutine hist(n, x, ncells, min_val, max_val, counts) bind(C) - integer(c_long), intent(in) :: n, ncells + integer(c_size_t), intent(in) :: n, ncells real(c_double), intent(in) :: x(n), min_val, max_val integer(c_long), intent(out) :: counts(ncells) !f2py intent(hide) :: n - integer(c_long) :: i, idx + integer(c_size_t) :: i, idx real(c_double) :: bin_width counts = 0_c_long @@ -181,11 +182,11 @@ subroutine hist(n, x, ncells, min_val, max_val, counts) bind(C) bin_width = 1._c_double ! prevent 0 division end if - do i=1, n + do i=1_c_size_t, n if (.NOT. isnan(x(i))) then - idx = int((x(i) - min_val) / bin_width, c_long) + 1 + idx = int((x(i) - min_val) / bin_width, c_size_t) + 1_c_size_t if (idx > ncells) then - idx = ncells + idx = int(ncells, c_size_t) end if counts(idx) = counts(idx) + 1 @@ -207,7 +208,7 @@ subroutine hist(n, x, ncells, min_val, max_val, counts) bind(C) ! counts : integer(8), counts in each bin ! -------------------------------------------------------------------- subroutine histogram(n, k, x, descriptor, counts) bind(C) - integer(c_long), intent(in) :: n, k + integer(c_size_t), intent(in) :: n, k real(c_double), intent(in) :: x(n) real(c_double), intent(out) :: descriptor(3) integer(c_long), intent(out) :: counts(k) diff --git a/src/skdh/io/_extensions/pyread.c b/src/skdh/io/_extensions/pyread.c index 6012f0bc..ffd99b77 100644 --- a/src/skdh/io/_extensions/pyread.c +++ b/src/skdh/io/_extensions/pyread.c @@ -10,6 +10,10 @@ #define NP_FROM_ANY(x) PyArray_FromAny(x, PyArray_DescrFromType(NPY_LONG), 1, 0, NPY_ARRAY_ENSUREARRAY | NPY_ARRAY_CARRAY_RO, NULL) +// Have to do this here so that we can use the Py_ssize_t type +extern void axivity_read_header(Py_ssize_t *, char[], AX_Info_t *, int *); + + void geneactiv_set_error_message(int ierr) { switch (ierr) diff --git a/src/skdh/io/_extensions/read_axivity.f95 b/src/skdh/io/_extensions/read_axivity.f95 index 099c74d1..a41d54db 100644 --- a/src/skdh/io/_extensions/read_axivity.f95 +++ b/src/skdh/io/_extensions/read_axivity.f95 @@ -97,13 +97,14 @@ subroutine axivity_close(info) bind(C, name="axivity_close") ! axivity_read_header : read the header (first 1024 bytes) plus a little bit for error checking ! ============================================================================================= subroutine axivity_read_header(flen, file, finfo, ierr) bind(C, name="axivity_read_header") - integer(c_long), intent(in) :: flen ! length of the file name string + integer(c_size_t), intent(in) :: flen ! length of the file name string character(kind=c_char), intent(in) :: file(flen) ! file name type(FileInfo_t), intent(inout) :: finfo ! file info storage structure integer(c_int), intent(inout) :: ierr ! error tracking/returning ! local type(metadata) :: hdr - integer(c_long) :: itmp, i + integer(c_size_t) :: i + integer(c_long) :: itmp character(960) :: annotation_block character(flen) :: file_ integer(c_int8_t) :: numAxesBps @@ -111,7 +112,7 @@ subroutine axivity_read_header(flen, file, finfo, ierr) bind(C, name="axivity_re integer(4) :: fstat_vals(13), fstat_err ! needed for converting c strings - do i=1_c_long, flen + do i=1_c_size_t, flen file_(i:i + 1) = file(i) end do diff --git a/src/skdh/io/_extensions/read_binary_imu.h b/src/skdh/io/_extensions/read_binary_imu.h index ab315d4e..85d1396c 100644 --- a/src/skdh/io/_extensions/read_binary_imu.h +++ b/src/skdh/io/_extensions/read_binary_imu.h @@ -74,7 +74,6 @@ typedef enum { AX_READ_E_BAD_LENGTH_ZERO_TIMESTAMPS = 7, } Read_Cwa_Error_t; -extern void axivity_read_header(long *, char[], AX_Info_t *, int *); extern void axivity_read_block(AX_Info_t *, long *, double *, double *, double *, int *); extern void adjust_timestamps(AX_Info_t *, double *, int *); extern void axivity_close(AX_Info_t *); diff --git a/src/skdh/io/_extensions/read_geneactiv.c b/src/skdh/io/_extensions/read_geneactiv.c index a40c249f..45459047 100644 --- a/src/skdh/io/_extensions/read_geneactiv.c +++ b/src/skdh/io/_extensions/read_geneactiv.c @@ -93,9 +93,9 @@ int get_timestamps(long *Nps, char time[40], GN_Info_t *info, GN_Data_t *data) data->ts[*Nps + j] = t0 + (double)j / info->fs; /* INDEXING */ - long mdays = MAX_DAYS; - long gns = GN_SAMPLES; - double block_t_delta = GN_SAMPLESf / info->fs; + // long mdays = MAX_DAYS; + // long gns = GN_SAMPLES; + // double block_t_delta = GN_SAMPLESf / info->fs; return GN_READ_E_NONE; } diff --git a/src/skdh/io/axivity.py b/src/skdh/io/axivity.py index 1300fa58..51b8bf68 100644 --- a/src/skdh/io/axivity.py +++ b/src/skdh/io/axivity.py @@ -29,6 +29,10 @@ class ReadCwa(BaseIO): of the key you do not want to trim. If provided, the tuple should be of the form (start_key, end_key). When provided, trim datetimes will be assumed to be in the same timezone as the data (ie naive if naive, or in the timezone provided). + trim_start_factor : int, optional + Factor of seconds to trim start time to. For example, if `trim_start_factor=15`, + then the start time will be trimmed to the next multiple of 15 seconds. Default + is None (no trimming). ext_error : {"warn", "raise", "skip"}, optional What to do if the file extension does not match the expected extension (.cwa). Default is "warn". "raise" raises a ValueError. "skip" skips the file @@ -48,14 +52,23 @@ class ReadCwa(BaseIO): {'accel': ..., 'time': ..., ...} """ - def __init__(self, *, trim_keys=None, ext_error="warn"): + def __init__(self, *, trim_keys=None, trim_start_factor=None, ext_error="warn"): super().__init__( # kwargs trim_keys=trim_keys, + trim_start_factor=trim_start_factor, ext_error=ext_error, ) - self.trim_keys = trim_keys + if trim_start_factor is not None: + if not 1 <= trim_start_factor <= 60: + raise ValueError("trim_start_factor must be between 1 and 60 seconds.") + + if trim_keys is None: + self.trim_keys = (None, None) + else: + self.trim_keys = trim_keys + self.trim_start_factor = trim_start_factor if ext_error.lower() in ["warn", "raise", "skip"]: self.ext_error = ext_error.lower() @@ -139,14 +152,14 @@ def predict(self, *, file, tz_name=None, **kwargs): if mag_axes is not None: # pragma: no cover :: don't have data to test this results[self._mag] = ascontiguousarray(imudata[:end, mag_axes]) - if self.trim_keys is not None: - results = self.trim_data( - *self.trim_keys, - tz_name, - kwargs, - **results, # contains the time array/argument - ) - + time, results = self.trim_data( + self.trim_start_factor, + *self.trim_keys, + tz_name, + kwargs, + **results, # contains the time array/argument + ) + results[self._time] = time results["fs"] = fs return results diff --git a/src/skdh/io/base.py b/src/skdh/io/base.py index 129db6c3..f3ddb206 100644 --- a/src/skdh/io/base.py +++ b/src/skdh/io/base.py @@ -9,7 +9,7 @@ import functools from warnings import warn -from numpy import nonzero +from numpy import nonzero, argwhere from pandas import to_datetime from skdh.base import BaseProcess @@ -133,7 +133,65 @@ def _to_timestamp(time_str, tz_name): return ts.timestamp() - def trim_data(self, start_key, end_key, tz_name, predict_kw, *, time, **data): + def trim_data(self, trim_start_factor, start_key, end_key, tz_name, predict_kw, *, time, **data): + """ + Trim data based on either a start factor or start/end keys + + Parameters + ---------- + trim_start_factor : int + Factor of seconds to trim start time to. For example, if `trim_start_factor=15`, + then the start time will be trimmed to the next multiple of 15 seconds. Computed + after trimming by time. + start_key : {None, str} + Start key for the provided trim start time. + end_key : {None, str} + End key for the provided trim end time. + tz_name : {None, str} + IANA time-zone name for the recording location. If not provided, + both the `time` array and provided trim times will be assumed to be naive + and in the same time-zone. + predict_kw : dict + Key-word arguments passed to the predict method. + """ + if start_key is not None or end_key is not None: + time, data = self.trim_data_time(start_key, end_key, tz_name, predict_kw, time=time, **data) + if trim_start_factor is not None: + time, data = self.trim_data_factor(trim_start_factor, time=time, **data) + + return time, data + + def trim_data_factor(self, trim_start_factor, *, time, **data): + """ + Trim raw data based on a factor of seconds. + + Parameters + ---------- + trim_start_factor : int + Factor of seconds to trim start time to. For example, if `trim_start_factor=15`, + then the start time will be trimmed to the next multiple of 15 seconds. + time : numpy.ndarray + Array of timestamps. + data : dict of numpy.ndarray + Dictionary of data arrays to trim, where each array is the same length as `time`. + """ + i1 = 0 + # convert the first N timestamps + for n in [1000, 5000, 10000, 50000]: + dt = to_datetime(time[:n], unit="s") + idx = argwhere((dt.second % trim_start_factor) == 0) + if idx.size > 0: + i1 = idx[0][0] + break + + self.logger.info(f"Trimming start time from {time[0]} to {time[i1]}, a factor of {trim_start_factor} seconds.") + + time = time[i1:] + res = {k: v[i1:] for k, v in data.items()} + + return time, res + + def trim_data_time(self, start_key, end_key, tz_name, predict_kw, *, time, **data): """ Trim raw data based on provided date-times @@ -153,11 +211,11 @@ def trim_data(self, start_key, end_key, tz_name, predict_kw, *, time, **data): trim_start = predict_kw.get(start_key) trim_end = predict_kw.get(end_key) - if trim_start is not None and trim_start is None: + if start_key is not None and trim_start is None: raise ValueError( f"`{start_key=}` was provided but not found in `predict` arguments" ) - if trim_end is not None and trim_end is None: + if end_key is not None and trim_end is None: raise ValueError( f"`{end_key=}` was provided but not found in `predict` arguments" ) @@ -178,7 +236,10 @@ def trim_data(self, start_key, end_key, tz_name, predict_kw, *, time, **data): i1 = nonzero(time <= max(time[0], ts_trim_start))[0][-1] i2 = nonzero(time >= min(time[-1], ts_trim_end))[0][0] - res = {self._time: time[i1:i2]} - res.update({k: v[i1:i2] for k, v in data.items()}) + self.logger.info(f"Trimming start time from {time[0]} to {time[i1]}") + self.logger.info(f"Trimming end time from {time[-1]} to {time[i2]}.") + + time = time[i1:i2] + res = {k: v[i1:i2] for k, v in data.items()} - return res + return time, res diff --git a/src/skdh/io/csv.py b/src/skdh/io/csv.py index f09b26f0..6f317428 100644 --- a/src/skdh/io/csv.py +++ b/src/skdh/io/csv.py @@ -48,6 +48,10 @@ class ReadCSV(BaseIO): of the key you do not want to trim. If provided, the tuple should be of the form (start_key, end_key). When provided, trim datetimes will be assumed to be in the same timezone as the data (ie naive if naive, or in the timezone provided). + trim_start_factor : int, optional + Factor of seconds to trim start time to. For example, if `trim_start_factor=15`, + then the start time will be trimmed to the next multiple of 15 seconds. Default + is None (no trimming). fill_gaps : bool, optional Fill any gaps in data streams. Default is True. If False and data gaps are detected, then the reading will raise a `ValueError`. @@ -123,6 +127,7 @@ def __init__( column_names, drop_duplicate_timestamps=False, trim_keys=None, + trim_start_factor=None, fill_gaps=True, fill_value=None, gaps_error="raise", @@ -138,6 +143,7 @@ def __init__( column_names=column_names, drop_duplicate_timestamps=drop_duplicate_timestamps, trim_keys=trim_keys, + trim_start_factor=trim_start_factor, fill_gaps=fill_gaps, fill_value=None, gaps_error=gaps_error, @@ -164,10 +170,18 @@ def __init__( if gaps_error.lower() not in ["raise", "warn", "ignore"]: raise ValueError("gaps_error must be one of `raise`, `warn`, or `ignore`.") + + if trim_start_factor is not None: + if not 1 <= trim_start_factor <= 60: + raise ValueError("trim_start_factor must be between 1 and 60 seconds.") self.time_col_name = time_col_name self.column_names = column_names - self.trim_keys = trim_keys + if trim_keys is None: + self.trim_keys = (None, None) + else: + self.trim_keys = trim_keys + self.trim_start_factor = trim_start_factor self.fill_gaps = fill_gaps self.fill_value = {} if fill_value is None else fill_value self.gaps_error = gaps_error @@ -364,9 +378,7 @@ def predict(self, *, file, tz_name=None, **kwargs): continue # trim the data - if self.trim_keys is not None: - data = self.trim_data(*self.trim_keys, tz_name, kwargs, time=time, **data) - time = data.pop(self._time) + time, data = self.trim_data(self.trim_start_factor, *self.trim_keys, tz_name, kwargs, time=time, **data) # now handle data gaps and second level timestamps, etc fs, time, results = self.handle_timestamp_inconsistency_np( diff --git a/src/skdh/io/geneactiv.py b/src/skdh/io/geneactiv.py index 5787b297..6fb13c70 100644 --- a/src/skdh/io/geneactiv.py +++ b/src/skdh/io/geneactiv.py @@ -38,14 +38,23 @@ class ReadBin(BaseIO): {'accel': ..., 'time': ...} """ - def __init__(self, trim_keys=None, ext_error="warn"): + def __init__(self, trim_keys=None, trim_start_factor=None, ext_error="warn"): super().__init__( # kwargs trim_keys=trim_keys, + trim_start_factor=trim_start_factor, ext_error=ext_error, ) - self.trim_keys = trim_keys + if trim_start_factor is not None: + if not 1 <= trim_start_factor <= 60: + raise ValueError("trim_start_factor must be between 1 and 60 seconds.") + + if trim_keys is None: + self.trim_keys = (None, None) + else: + self.trim_keys = trim_keys + self.trim_start_factor = trim_start_factor if ext_error.lower() in ["warn", "raise", "skip"]: self.ext_error = ext_error.lower() @@ -97,30 +106,21 @@ def predict(self, *, file, tz_name=None, **kwargs): n_max, fs, acc, time, light, temp = read_geneactiv(str(file)) # trim the data if necessary - if self.trim_keys is not None: - results = self.trim_data( - *self.trim_keys, - tz_name, - kwargs, - time=handle_naive_timestamps( - time[:n_max], is_local=True, tz_name=tz_name - ), - **{ - self._acc: acc[:n_max, :], - self._temp: temp[:n_max], - "light": light[:n_max], - }, - ) - results["fs"] = fs - else: - results = { - self._time: handle_naive_timestamps( - time[:n_max], is_local=True, tz_name=tz_name - ), + time, results = self.trim_data( + self.trim_start_factor, + *self.trim_keys, + tz_name, + kwargs, + time=handle_naive_timestamps( + time[:n_max], is_local=True, tz_name=tz_name + ), + **{ self._acc: acc[:n_max, :], self._temp: temp[:n_max], "light": light[:n_max], - "fs": fs, - } + }, + ) + results["fs"] = fs + results[self._time] = time return results diff --git a/src/skdh/utility/_extensions/meson.build b/src/skdh/utility/_extensions/meson.build index bf473c5c..5ab9fecb 100644 --- a/src/skdh/utility/_extensions/meson.build +++ b/src/skdh/utility/_extensions/meson.build @@ -15,7 +15,12 @@ movstat_lib = static_library( 'stack.c', 'moving_extrema.c', ], - c_args: numpy_nodepr_api, + c_args: [ + numpy_nodepr_api, + '-fno-unsafe-math-optimizations', + '-frounding-math', + '-fsignaling-nans', + ], include_directories: [inc_np], ) diff --git a/src/skdh/utility/_extensions/moving_moments.f95 b/src/skdh/utility/_extensions/moving_moments.f95 index 23fdc813..93db5283 100644 --- a/src/skdh/utility/_extensions/moving_moments.f95 +++ b/src/skdh/utility/_extensions/moving_moments.f95 @@ -3,6 +3,42 @@ ! Copyright (c) 2021. Pfizer Inc. All rights reserved. +subroutine mov_moments_1_full(n, x, wlen, skip, req_n, mean) bind(C, name="mov_moments_1_full") + use, intrinsic :: iso_c_binding + use, intrinsic :: ieee_arithmetic, only: IEEE_Value, IEEE_QUIET_NAN + implicit none + integer(c_long), intent(in) :: n, wlen, skip, req_n + real(c_double), intent(in) :: x(n) + real(c_double), intent(out) :: mean((n - 1)/skip + 1) + ! local + integer(c_long) :: i, j + real(c_double) :: m1(n) + + m1(1) = x(1) + + do i=2, n + m1(i) = m1(i-1) + x(i) + end do + + j = 2_c_long + mean(1) = m1(wlen) + + do i=skip, n-1, skip + if ((i + wlen) < n) then + mean(j) = m1(i+wlen) - m1(i) + else + if ((n - i) >= req_n) then + mean(j) = (m1(n) - m1(i)) * wlen / real(n - i, c_double) + else + mean(j) = IEEE_Value(mean(j), IEEE_QUIET_NAN) + end if + end if + j = j + 1 + end do + + mean = mean / wlen +end subroutine + subroutine mov_moments_1(n, x, wlen, skip, mean) bind(C, name="mov_moments_1") use, intrinsic :: iso_c_binding implicit none diff --git a/src/skdh/utility/_extensions/moving_statistics.c b/src/skdh/utility/_extensions/moving_statistics.c index 1a72e958..b64c81f6 100644 --- a/src/skdh/utility/_extensions/moving_statistics.c +++ b/src/skdh/utility/_extensions/moving_statistics.c @@ -13,6 +13,7 @@ /* moving moments */ extern void mov_moments_1(long *, double *, long *, long *, double *); +extern void mov_moments_1_full(long *, double *, long *, long *, long *, double *); extern void moving_moments_1(long *, double *, long *, long *, double *); extern void mov_moments_2(long *, double *, long *, long *, double *, double *); extern void moving_moments_2(long *, double *, long *, long *, double *, double *); @@ -24,10 +25,9 @@ extern void fmoving_median(long *, double *, long *, long *, double *); PyObject * moving_mean(PyObject *NPY_UNUSED(self), PyObject *args){ PyObject *x_; - long wlen, skip; - int trim; + long wlen, skip, req_npts; - if (!PyArg_ParseTuple(args, "Ollp:moving_mean", &x_, &wlen, &skip, &trim)) + if (!PyArg_ParseTuple(args, "Olll:moving_mean", &x_, &wlen, &skip, &req_npts)) return NULL; PyArrayObject *data = (PyArrayObject *)PyArray_FromAny( @@ -45,7 +45,6 @@ PyObject * moving_mean(PyObject *NPY_UNUSED(self), PyObject *args){ int ndim = PyArray_NDIM(data); const npy_intp *ddims = PyArray_DIMS(data); long npts = ddims[ndim - 1]; - long trim_pts = (npts - wlen) / skip + 1; npy_intp *rdims = (npy_intp *)malloc(ndim * sizeof(npy_intp)); if (!rdims) { @@ -58,20 +57,15 @@ PyObject * moving_mean(PyObject *NPY_UNUSED(self), PyObject *args){ rdims[i] = ddims[i]; } // dimension of the roll - if (trim) - { - rdims[ndim - 1] = trim_pts; - } else { - rdims[ndim - 1] = (npts - 1) / skip + 1; - } + rdims[ndim - 1] = (npts - 1) / skip + 1; PyArrayObject *rmean = (PyArrayObject *)PyArray_EMPTY(ndim, rdims, NPY_DOUBLE, 0); - free(rdims); if (!rmean) { Py_XDECREF(data); Py_XDECREF(rmean); + free(rdims); return NULL; } @@ -84,16 +78,13 @@ PyObject * moving_mean(PyObject *NPY_UNUSED(self), PyObject *args){ for (int i = 0; i < nrepeats; ++i) { - for (int j = trim_pts; j < res_stride; ++j) - { - rmean_ptr[j] = NPY_NAN; - } - mov_moments_1(&npts, dptr, &wlen, &skip, rmean_ptr); + mov_moments_1_full(&npts, dptr, &wlen, &skip, &req_npts, rmean_ptr); dptr += npts; // increment by number of points in last dimension rmean_ptr += res_stride; } Py_XDECREF(data); + free(rdims); return (PyObject *)rmean; } diff --git a/src/skdh/utility/internal.py b/src/skdh/utility/internal.py index 22e75315..43fada66 100644 --- a/src/skdh/utility/internal.py +++ b/src/skdh/utility/internal.py @@ -27,6 +27,7 @@ concatenate, roll, full, + minimum, ) from scipy.signal import cheby1, sosfiltfilt @@ -169,8 +170,8 @@ def get_day_index_intersection(starts, stops, for_inclusion, day_start, day_stop if fs is None: raise ValueError("Sampling frequency `fs` must be provided if `ends_round` is provided.") factor = int(ends_round * fs) - valid_starts = around((valid_starts - day_start) / factor, 0) * factor + day_start - valid_stops = around((valid_stops - day_start) / factor, 0) * factor + day_start + valid_starts = minimum(around((valid_starts - day_start) / factor, 0) * factor + day_start, day_stop) + valid_stops = minimum(around((valid_stops - day_start) / factor, 0) * factor + day_start, day_stop) # make sure they stay as integers valid_starts = valid_starts.astype(int_) diff --git a/src/skdh/utility/math.py b/src/skdh/utility/math.py index 12027908..dc6c00cb 100644 --- a/src/skdh/utility/math.py +++ b/src/skdh/utility/math.py @@ -42,7 +42,7 @@ ] -def moving_mean(a, w_len, skip, trim=True, axis=-1): +def moving_mean(a, w_len, skip, req_points=1.0, trim=True, axis=-1): r""" Compute the moving mean. @@ -54,6 +54,9 @@ def moving_mean(a, w_len, skip, trim=True, axis=-1): Window length in number of samples. skip : int Window start location skip in number of samples. + req_points : float, optional + Minimum fraction of points required to compute the mean. Default is 1.0. + If a float, is a proportion of the window length. If an int, is the number of points. trim : bool, optional Trim the ends of the result, where a value cannot be calculated. If False, these values will be set to NaN. Default is True. @@ -132,6 +135,11 @@ def moving_mean(a, w_len, skip, trim=True, axis=-1): """ if w_len <= 0 or skip <= 0: raise ValueError("`wlen` and `skip` cannot be less than or equal to 0.") + + if isinstance(req_points, float): + req_points = int(round(w_len * req_points)) + if req_points < 0: + raise ValueError("`req_points` cannot be less than 0.") # move computation axis to end x = moveaxis(a, axis, -1) @@ -140,7 +148,12 @@ def moving_mean(a, w_len, skip, trim=True, axis=-1): if w_len > x.shape[-1]: raise ValueError("Window length is larger than the computation axis.") - rmean = _extensions.moving_mean(x, w_len, skip, trim) + rmean = _extensions.moving_mean(x, w_len, skip, req_points) + + # handle trim + if trim: + n = (x.shape[-1] - w_len) // skip + 1 + rmean = rmean[..., :n] # move computation axis back to original place and return return moveaxis(rmean, -1, axis) diff --git a/tests/utility/test_math.py b/tests/utility/test_math.py index 112cf5c8..82c99d58 100644 --- a/tests/utility/test_math.py +++ b/tests/utility/test_math.py @@ -160,6 +160,29 @@ class TestMovingMean(BaseMovingStatsTester): truth_function = staticmethod(mean) truth_kw = {} + @pytest.mark.parametrize("skip", (1, 2, 150, 300)) + @pytest.mark.parametrize("req_points", (1, 0.5, 250)) + def test_full(self, skip, req_points, np_rng): + wlen = 250 # constant + x = np_rng.random(2000) + + if isinstance(req_points, float): + t_req_points = int(round(req_points * wlen)) + else: + t_req_points = req_points + truth = full((x.size - 1) // skip + 1, nan) + j = 0 + for i in range(0, x.size, skip): + if (x.size - i) >= t_req_points: + truth[j] = mean(x[i:i + wlen]) + else: + truth[j] = nan + j += 1 + + pred = self.function(x, wlen, skip, trim=False, req_points=req_points) + + assert allclose(pred, truth, equal_nan=True) + class TestMovingSD(BaseMovingStatsTester): function = staticmethod(moving_sd)