R/mims_unit.R

In MIMSunit: Algorithm to Compute Monitor Independent Movement Summary Unit (MIMS-Unit)

#' Compute Monitor Independent Motion Summary unit (MIMS-unit)
#'
#' \code{mims_unit} computes the Monitor Independent Motion Summary unit for the
#' input multi-channel accelerometer signal. The input signal can be from
#' devices of any sampling rate and dynamic range. Please refer to the
#' manuscript for detailed description of the algorithm. Please refer to
#' functions for the intermediate steps: \code{\link{extrapolate}} for
#' extrapolation, \code{\link{iir}} for filtering,
#' \code{\link{aggregate_for_mims}} for aggregation.
#'
#' @note This function is a wrapper function for the low-level
#'   \code{\link{custom_mims_unit}} function. It has set internal parameters as
#'   described in the manuscript. If you want to run customized algorithm for
#'   MIMSunit or if you want to develop better algorithms based on MIMS-unit
#'   algorithm, please use function \code{\link{custom_mims_unit}} where all
#'   parameters are tunable.
#'
#'   \code{before_df} and \code{after_df} are often set when the accelerometer
#'   data are divided into files of smaller chunk.
#'
#'   Please make sure the input data do not contain duplicated timestamps. See
#'    more information about this \href{https://github.com/mHealthGroup/MIMSunit/issues/32}{issue}.
#'    Otherwise the computation will stop.
#'
#' @section How is it used in MIMS-unit algorithm?: This is the main entry of
#'   MIMS-unit algorithm.
#'
#' @param epoch string. Any format that is acceptable by argument \code{breaks}
#'   in method \code{\link[base]{cut.POSIXt}}.For example, "1 sec", "1 min", "5
#'   sec", "10 min". Default is "5 sec".
#' @param dynamic_range numerical vector. The dynamic ranges of the input
#'   signal. Should be a 2-element numerical vector. \code{c(low, high)}, where
#'   \code{low} is the negative max value the device can reach and \code{high}
#'   is the positive max value the device can reach.
#' @param output_mims_per_axis logical. If it is TRUE, the output MIMS-unit
#'   dataframe will have MIMS-unit values for each axis from the third column.
#'   Default is FALSE.
#' @param use_gui_progress logical. If it is TRUE, show GUI progress bar on
#'   windows platform. Default is FALSE.
#' @param use_snapshot_to_check logical. If TRUE, the function will use the first
#'   100 rows or 10% (whichever is smaller) to check timestamp duplications. Otherwise,
#'   the algorithm will use all data to check timestamp duplications. Default is FALSE.
#' @return dataframe. The MIMS-unit dataframe. The first column is the start
#'   time of each epoch in POSIXct format. The second column is the MIMS-unit
#'   value for the input signal. If \code{output_mims_per_axis} is TRUE, the
#'   third column and then are the MIMS-unit values for each axis of the input
#'   signal.
#'
#' @family Top level API functions
#' @name mims_unit


#' @rdname mims_unit
#' @param df dataframe. Input multi-channel accelerometer signal. The
#' first column should be the time component. The accelerometer data values
#'  (typically starting from the second column) should be in \eqn{g} (per \eqn{9.81m/s^2}) unit.
#' @param before_df dataframe. The multi-channel accelerometer signal comes
#'   before the input signal to be prepended to the input signal during
#'   computation. This is used to eliminate the edge effect during extrapolation
#'   and filtering. If it is \code{NULL}, algorithm will run directly on the
#'   input signal. Default is NULL.
#' @param after_df dataframe. The multi-channel accelerometer signal comes after
#'   the input signal to be append to the input signal. This is used to
#'   eliminate the edge effect during extrapolation and filtering. If it is
#'   \code{NULL}, algorithm will run directly on the input signal. Default is
#'   NULL.
#' @param st character or POSIXct timestamp. An optional start time you can set to
#'   force the epochs generated by referencing this start time. If it is NULL, the
#'   function will use the first timestamp in the timestamp column as start time to
#'   generate epochs. This is useful when you are processing a stream of data and
#'   want to use a common start time for segmenting data. Default is NULL.
#' @export
#' @examples
#'   # Use sample data for testing
#'   df = sample_raw_accel_data
#'
#'   # compute mims unit values and output axial values
#'   output = mims_unit(df, epoch = '1 sec', dynamic_range=c(-8, 8), output_mims_per_axis=TRUE)
#'   head(output)
mims_unit <-
  function(df,
           before_df = NULL,
           after_df = NULL,
           epoch = "5 sec",
           dynamic_range,
           output_mims_per_axis = FALSE,
           use_gui_progress = FALSE,
           st = NULL,
           use_snapshot_to_check=FALSE) {
    mims_df <- custom_mims_unit(
      df = df,
      epoch = epoch,
      dynamic_range = dynamic_range,
      noise_level = 0.03,
      k = 0.05,
      spar = 0.6,
      filter_type = "butter",
      cutoffs = c(0.2, 5),
      axes = c(2, 3, 4),
      use_extrapolation = TRUE,
      use_filtering = TRUE,
      combination = "sum",
      allow_truncation = TRUE,
      output_mims_per_axis = output_mims_per_axis,
      output_orientation_estimation = FALSE,
      before_df = before_df,
      after_df = after_df,
      use_gui_progress = use_gui_progress,
      st = st,
      use_snapshot_to_check = use_snapshot_to_check
    )
    return(mims_df)
  }

#' Estimates sensor orientation
#'
#' \code{sensor_orientations} estimates the orientation angles for the input
#' multi-channel accelerometer signal. The input signal can be from devices of
#' any sampling rate and dynamic range. Please refer to function
#' \code{\link{compute_orientation}} for the implementation of the estimation
#' algorithm.
#'
#' @note This function interpolates and extrapolates the signal before
#'   estimating the orientation angles.
#'
#'   \code{before_df} and \code{after_df} are often set when the accelerometer
#'   data are divided into files of smaller chunk.
#'
#' @section How is it used in MIMS-unit algorithm?: This is not included in the
#'   official MIMS-unit algorithm nor the manuscript, but we found it is useful
#'   to know the sensor orientations in addition to the summary of movement.
#'
#' @param df dataframe. Input multi-channel accelerometer signal.
#' @param before_df dataframe. The multi-channel accelerometer signal comes
#'   before the input signal to be prepended to the input signal during
#'   computation. This is used to eliminate the edge effect during extrapolation
#'   and filtering. If it is \code{NULL}, algorithm will run directly on the
#'   input signal. Default is NULL.
#' @param after_df dataframe. The multi-channel accelerometer signal comes after
#'   the input signal to be append to the input signal. This is used to
#'   eliminate the edge effect during extrapolation and filtering. If it is
#'   \code{NULL}, algorithm will run directly on the input signal. Default is
#'   NULL.
#' @param epoch string. Any format that is acceptable by argument \code{breaks}
#'   in method \code{\link[base]{cut.POSIXt}}.For example, "1 sec", "1 min", "5
#'   sec", "10 min". Default is "5 sec".
#' @param dynamic_range numerical vector. The dynamic ranges of the input
#'   signal. Should be a 2-element numerical vector. \code{c(low, high)}, where
#'   \code{low} is the negative max value the device can reach and \code{high}
#'   is the positive max value the device can reach.
#' @param st character or POSIXct timestamp. An optional start time you can set to
#'   force the epochs generated by referencing this start time. If it is NULL, the
#'   function will use the first timestamp in the timestamp column as start time to
#'   generate epochs. This is useful when you are processing a stream of data and
#'   want to use a common start time for segmenting data. Default is NULL.
#' @return dataframe. The orientation dataframe. The first column is the start
#'   time of each epoch in POSIXct format. The second to fourth columns are the
#'   orientation angles.
#'
#' @family Top level API functions
#' @export
#' @examples
#'   # Use sample data for testing
#'   df = sample_raw_accel_data
#'
#'   # compute sensor orientation angles
#'   sensor_orientations(df, epoch = '2 sec', dynamic_range=c(-8, 8))
#'
#'   # compute sensor orientation angles with different epoch length
#'   output = sensor_orientations(df, epoch = '1 sec', dynamic_range=c(-8, 8))
#'   head(output)
sensor_orientations <-
  function(df,
           before_df = NULL,
           after_df = NULL,
           epoch = "5 sec",
           dynamic_range,
           st = NULL) {
    ori_df <- custom_mims_unit(
      df = df,
      epoch = epoch,
      dynamic_range = dynamic_range,
      output_orientation_estimation = TRUE,
      epoch_for_orientation_estimation = epoch,
      before_df = before_df,
      after_df = after_df,
      st = st
    )[[2]]
    return(ori_df)
  }


#' Compute both MIMS-unit and sensor orientations with custom settings
#'
#' \code{custom_mims_unit} computes the Monitor Independent Motion Summary unit
#' and estimates the sensor orientations for the input multi-channel
#' accelerometer signal with custom settings. The input signal can be from
#' devices of any sampling rate and dynamic range. Please refer to the
#' manuscript for detailed description of the algorithm. Please refer to
#' functions for the intermediate steps: \code{\link{extrapolate}} for
#' extrapolation, \code{\link{iir}} for filtering,
#' \code{\link{aggregate_for_mims}} and \code{\link{aggregate_for_orientation}}
#' for aggregation.
#'
#' @note This function allows you to run customized algorithm for MIMSunit and
#'   sensor orientations.
#'
#'   \code{before_df} and \code{after_df} are often set when the accelerometer
#'   data are divided into files of smaller chunk.
#'
#' @section How is it used in MIMS-unit algorithm?: This is the low-level entry
#'   of MIMS-unit and orientation estimation algorithm. \code{\link{mims_unit}}
#'   calls this function internally.
#'
#' @param df dataframe. Input multi-channel accelerometer signal.
#' @param epoch string. Any format that is acceptable by argument \code{breaks}
#'   in method \code{\link[base]{cut.POSIXt}}.For example, "1 sec", "1 min", "5
#'   sec", "10 min". Default is "5 sec".
#' @param dynamic_range numerical vector. The dynamic ranges of the input
#'   signal. Should be a 2-element numerical vector. \code{c(low, high)}, where
#'   \code{low} is the negative max value the device can reach and \code{high}
#'   is the positive max value the device can reach.
#' @param noise_level number. The tolerable noise level in \eqn{g} unit, should
#'   be between 0 and 1. Default is 0.03, which applies to most devices.
#' @param k number. Duration of neighborhood to be used in local spline
#'   regression for each side, in seconds. Default is 0.05, as optimized by
#'   MIMS-unit algorithm.
#' @param spar number. Between 0 and 1, to control how smooth we want to fit
#'   local spline regression, 0 is linear and 1 matches all local points.
#'   Default is 0.6, as optimized by MIMS-unit algorithm.
#' @param filter_type string. The type of filter to be applied. Could be
#'   'butter' for butterworth bandpass filter, 'ellip' for elliptic bandpass
#'   filter or 'bessel' for bessel lowpass filter + average removal highpass
#'   filter. Default is "butter".
#' @param cutoffs numerical vector. Cut off frequencies to be used in filtering.
#'   If \code{filter_type} is "bessel", the cut off frequency for lowpass filter
#'   would be multiplied by 2 when being used. Default is 0.2Hz and 5Hz.
#' @param axes numerical vector. Indices of columns that specifies the axis
#'   values of the input signal. Default is \code{c(2,3,4)}.
#' @param use_extrapolation logical. If it is TRUE, the function will apply
#'   extrapolation algorithm to the input signal, otherwise it will skip
#'   extrapolation but only linearly interpolate the signal to 100Hz. Default is
#'   TRUE.
#' @param use_filtering logical. If it is TRUE, the function will apply bandpass
#'   filtering to the input signal, otherwise it will skip the filtering.
#'   Default is TRUE.
#' @param combination string. Method to combine MIMS-unit values for each axis.
#'   Could be "sum" for \code{\link{sum_up}} or "vm" for
#'   \code{\link{vector_magnitude}}.
#' @param allow_truncation logical. If it is TRUE, the algorithm will truncate
#'   very small MIMS-unit values to zero. Default is TRUE.
#' @param output_mims_per_axis logical. If it is TRUE, the output MIMS-unit
#'   dataframe will have MIMS-unit values for each axis from the third column.
#'   Default is FALSE.
#' @param output_orientation_estimation logical. If it is TRUE, the function
#'   will also estimate sensor orientations over each epoch. And the output will
#'   be a list, with the first element being the MIMS-unit dataframe, and the
#'   second element being the sensor orientation dataframe. Default is FALSE.
#' @param epoch_for_orientation_estimation string. string. Any format that is
#'   acceptable by argument \code{breaks} in method
#'   \code{\link[base]{cut.POSIXt}}.For example, "1 sec", "1 min", "5 sec", "10
#'   min". Default is "5 sec". It is independent from \code{epoch} for
#'   MIMS-unit.
#' @param before_df dataframe. The multi-channel accelerometer signal comes
#'   before the input signal to be prepended to the input signal during
#'   computation. This is used to eliminate the edge effect during extrapolation
#'   and filtering. If it is \code{NULL}, algorithm will run directly on the
#'   input signal. Default is NULL.
#' @param after_df dataframe. The multi-channel accelerometer signal comes after
#'   the input signal to be append to the input signal. This is used to
#'   eliminate the edge effect during extrapolation and filtering. If it is
#'   \code{NULL}, algorithm will run directly on the input signal. Default is
#'   NULL.
#' @param use_gui_progress logical. If it is TRUE, show GUI progress bar on
#'   windows platform. Default is FALSE.
#' @param st character or POSIXct timestamp. An optional start time you can set to
#'   force the epochs generated by referencing this start time. If it is NULL, the
#'   function will use the first timestamp in the timestamp column as start time to
#'   generate epochs. This is useful when you are processing a stream of data and
#'   want to use a common start time for segmenting data. Default is NULL.
#' @param use_snapshot_to_check logical. If TRUE, the function will use the first
#'   100 rows or 10% (whichever is smaller) to check timestamp duplications. Otherwise,
#'   the algorithm will use all data to check timestamp duplications. Default is FALSE.
#' @return dataframe or list. If \code{output_orientation_estimation} is TRUE,
#'   the output will be a list, otherwise the output will be the MIMS-unit
#'   dataframe.
#'
#'   The first element will be the MIMS-unit dataframe, in which the first
#'   column is the start time of each epoch in POSIXct format, and the second
#'   column is the MIMS-unit value for the input signal, and the third column
#'   and on are the MIMS-unit values for each axis of the input signal if
#'   \code{output_mims_per_axis} is TRUE.
#'
#'   The second element will be the orientation dataframe, in which the first
#'   column is the start time of each epoch in POSIXct format, and the second to
#'   fourth column is the estimated orientations for the input signal.
#'
#' @family Top level API functions
#' @export
#' @examples
#'   # Use sample data for testing
#'   df = sample_raw_accel_data
#'
#'   # compute mims unit values with custom parameter
#'   output = custom_mims_unit(df, epoch = '1 sec', dynamic_range=c(-8, 8), spar=0.7)
#'   head(output)
custom_mims_unit <-
  function(df,
           epoch = "5 sec",
           dynamic_range,
           noise_level = 0.03,
           k = 0.05,
           spar = 0.6,
           filter_type = "butter",
           cutoffs = c(0.2, 5),
           axes = c(2, 3, 4),
           use_extrapolation = TRUE,
           use_filtering = TRUE,
           combination = "sum",
           allow_truncation = TRUE,
           output_mims_per_axis = FALSE,
           output_orientation_estimation = FALSE,
           epoch_for_orientation_estimation = NULL,
           before_df = NULL,
           after_df = NULL,
           use_gui_progress = FALSE,
           st = NULL,
           use_snapshot_to_check=FALSE) {


    if (inherits(df, "tbl_df")) {
      df = as.data.frame(df)
    }
    first_col = df[[1]]
    if (is.unsorted(first_col)) {
      df = df[ order(first_col), ]
    }
    rm(first_col)

    # check timestamp duplication after the timestamp column is sorted
    if (use_snapshot_to_check) {
      ts_has_duplication = any(diff(df[1:pmax(100,round(nrow(df[,1])*0.1)),1]) == 0)
    } else {
      ts_has_duplication = any(diff(df[,1]) == 0)
    }
    if (ts_has_duplication) {
      stop("Input data contains duplicated timestamps!")
    }

    # save the start and stop time of original df
    if (shiny::isRunning()) {
      pb = shiny::Progress$new()
      pb$set(message = "Computing MIMS-unit values", value=0)

    }
    else if (use_gui_progress & .Platform$OS.type == 'windows') {
      ProgressBar = utils::winProgressBar
      setProgressBar = utils::setWinProgressBar
      pb = ProgressBar(min = 0,
                       max = 1,
                       title = "Computing MIMS-unit values",
                       label = 'Starting...')
    }
    else {
      ProgressBar = utils::txtProgressBar
      setProgressBar = utils::setTxtProgressBar
      pb = ProgressBar(min = 0,
                       max = 1,
                       title = "Computing MIMS-unit values",
                       label = 'Starting...')
    }



    start_time <- lubridate::floor_date(df[[1]][1], unit = "seconds")
    stop_time <-
      lubridate::floor_date(df[[1]][nrow(df)], unit = "seconds")

    # concatenate with before and after df
    if (is.data.frame(before_df)) {
      if (shiny::isRunning()) {
        pb$inc(0.2, detail='Concatenating before_df...')
      } else {
        setProgressBar(pb, 0.2, label = 'Concatenating before_df...')
      }
      df <- rbind(before_df, df)
    }
    if (is.data.frame(after_df)) {
      if (shiny::isRunning()) {
        pb$inc(0.3, detail='Concatenating after_df...')
      } else {
        setProgressBar(pb, 0.3, label = 'Concatenating after_df...')
      }
      df <- rbind(df, after_df)
    }

    # apply extrapolation algorithm
    if (use_extrapolation) {
      stopifnot(length(dynamic_range) == 2)
      if (shiny::isRunning()) {
        pb$inc(0.4, detail='Extrapolating signal...')
      } else {
        setProgressBar(pb, 0.4, label = 'Extrapolating signal...')
      }
      resampled_data <-
        extrapolate(df, dynamic_range, noise_level, k, spar)
    } else {
      if (shiny::isRunning()) {
        pb$inc(0.4, detail='Interpolating signal...')
      } else {
        setProgressBar(pb, 0.4, label = 'Interpolating signal...')
      }
      resampled_data <-
        interpolate_signal(df, sr = 100, method = "linear")
    }
    rm(df)
    # removed the extrapolated_data
    sr <- sampling_rate(resampled_data)


    # store -0.01 values separately
    row_abnormal <- rep(FALSE, nrow(resampled_data))
    for (i in 2:ncol(resampled_data))
    {
      row_abnormal <- row_abnormal | resampled_data[[i]] < -150
    }

    abnormal_data <- resampled_data[row_abnormal, ]
    normal_data <- resampled_data[!row_abnormal, ]

    # Compute orientations
    if (output_orientation_estimation) {
      if (is.null(epoch_for_orientation_estimation)) {
        epoch_for_orientation_estimation <- epoch
      }
      orientation_data <-
        aggregate_for_orientation(resampled_data,
                                  epoch = epoch_for_orientation_estimation,
                                  st = st
        )
    } else {
      orientation_data <- NULL
    }
    rm(resampled_data); gc();

    # Apply filter cascade
    if (use_filtering) {
      if (shiny::isRunning()) {
        pb$inc(0.6, detail='Filtering signal...')
      } else {
        setProgressBar(pb, 0.6, label = 'Filtering signal...')
      }
      if (filter_type == "butter") {
        filtered_data <-
          iir(
            normal_data,
            sr = sr,
            cutoff_freq = cutoffs,
            order = 4,
            type = "pass",
            filter_type = "butter"
          )
      } else if (filter_type == "ellip") {
        filtered_data <-
          iir(
            normal_data,
            sr = sr,
            cutoff_freq = cutoffs,
            order = 4,
            type = "pass",
            filter_type = "ellip"
          )
      }
    } else {
      filtered_data <- normal_data
    }
    rm(normal_data); gc()

    # sort by timestamp
    colnames(abnormal_data)[2:4] <- colnames(filtered_data)[2:4]
    filtered_data <- rbind(filtered_data, abnormal_data)
    rm(abnormal_data)
    filtered_data <-
      filtered_data[order(filtered_data[[1]]), ]


    # Compute the AUC
    if (shiny::isRunning()) {
      pb$inc(0.8, detail='Computing AUC...')
    } else {
      setProgressBar(pb, 0.8, label = 'Computing AUC...')
    }
    integrated_data <-
      aggregate_for_mims(
        filtered_data,
        epoch = epoch,
        method = "trapz",
        rectify = TRUE,
        st = st
      )
    rm(filtered_data);

    if (allow_truncation) {
      if (shiny::isRunning()) {
        pb$inc(0.9, detail='Truncating signal...')
      } else {
        setProgressBar(pb, 0.9, label = 'Truncating signal...')
      }
      truncate_indices <-
        integrated_data[, 2:ncol(integrated_data)] > 0 &
          (integrated_data[, 2:ncol(integrated_data)] <=
            (1e-04 * parse_epoch_string(epoch, sr)))
      truncate_indices <- data.frame(truncate_indices)
      integrated_data[, 2:ncol(integrated_data)] <-
        sapply(1:(ncol(integrated_data) - 1), function(n) {
          integrated_data[truncate_indices[, n], n + 1] <- 0
          return(integrated_data[, n + 1])
        })
    }

    # Compute vector magnitude
    if (shiny::isRunning()) {
      pb$inc(1, detail='Summing up axial values...')
      pb$close()
    } else {
      setProgressBar(pb, 1, label = 'Summing up axial values...')
      close(pb)
    }
    row_abnormal <- rep(FALSE, nrow(integrated_data))
    for (i in 2:ncol(integrated_data))
    {
      row_abnormal <- row_abnormal | integrated_data[[i]] < 0
    }
    if (combination == "vm") {
      mims_data <-
        vector_magnitude(integrated_data, axes = axes)
    } else if (combination == "sum") {
      mims_data <- sum_up(integrated_data, axes = axes)
    } else {
      mims_data <- sum_up(integrated_data, axes = axes)
    }

    if (output_mims_per_axis) {
      mims_data <-
        cbind(mims_data, integrated_data[, 2:ncol(integrated_data)])
      colnames(mims_data)[2:ncol(mims_data)] <-
        c("MIMS_UNIT", "MIMS_UNIT_X", "MIMS_UNIT_Y", "MIMS_UNIT_Z")
      mims_data[row_abnormal, c(2, 3, 4, 5)] <- -0.01
    } else {
      colnames(mims_data)[2] <- "MIMS_UNIT"
      mims_data[row_abnormal, 2] <- -0.01
    }
    rm(row_abnormal)
    rm(integrated_data)


    # only keep data between start and end time
    keep_mask <-
      mims_data[[1]] >= start_time & mims_data[[1]] < stop_time
    mims_data <- mims_data[keep_mask, ]


    if (output_orientation_estimation) {
      return(list(mims = mims_data, orientation = orientation_data))
    } else {
      return(mims_data)
    }
  }

#' @rdname mims_unit
#'
#' @param files character vector. A list of csv filepaths for raw accelerometer
#'   data organized in order to be processed. The data should be consecutive in
#'   timestamps. A typical case is a set of hourly or daily files for
#'   continuous accelerometer sampling. For a single file, please wrap the filepath
#'   in a vector `c(filepath)`.
#' @param file_type character. "mhealth" or "actigraph". The type of the csv files
#' that store the raw accelerometer data.
#' @param ... additional parameters passed to the import function when reading
#'   in the data from the files.
#' @export
#' @examples
#'   # Use sample mhealth file for testing
#'   filepaths = c(
#'     system.file('extdata', 'mhealth.csv', package='MIMSunit')
#'   )
#'
#'   # Test with multiple files
#'   output = mims_unit_from_files(filepaths, epoch = "1 sec", dynamic_range = c(-8, 8))
#'   head(output)
mims_unit_from_files <-
  function(files,
           epoch = "5 sec",
           dynamic_range,
           output_mims_per_axis = FALSE,
           use_gui_progress = FALSE,
           use_snapshot_to_check = FALSE,
           file_type = "mhealth", ...) {
    num_of_files <- length(files)
    dots = list(...)
    df <- NULL
    results <- list()

    if (file_type == "mhealth") {
      import_fun <- import_mhealth_csv_chunked
      header = TRUE
    } else if (file_type == "actigraph") {
      import_fun <- function(x, chunk_samples) import_actigraph_csv_chunked(x, chunk_samples = chunk_samples, ...)
      header = dots$header
    } else {
      stop('Only "mhealth" or "actigraph" file types are supported')
    }

    meta = .get_meta(files[1], file_type = file_type, header = header)
    sr = meta$sr
    st = meta$st
    num_samples_epoch = parse_epoch_string(epoch, sr)
    num_per_load = max(num_samples_epoch, 180000)
    last_epoch_st = NULL
    last_chunk = NULL
    j = 1
    if (shiny::isRunning()) {
      file_pb = shiny::Progress$new()
    }
    for (i in 1:num_of_files) {
      if (shiny::isRunning()) {
        file_pb$set(value=i/num_of_files, message=paste("Compute MIMS-unit values for", files[i]))
      }
      funcs = import_fun(files[i], chunk_samples = num_per_load)
      next_chunk = funcs[[1]]
      close_con = funcs[[2]]

      if (shiny::isRunning()) {
        chunk_pb = shiny::Progress$new()
        chunk_pb$set(value=0, message="Compute MIMS-unit values for chunks")
      }
      repeat {
        chunk = next_chunk()
        if (nrow(chunk) > 0) {
          if (shiny::isRunning()) {
            chunk_pb$inc(amount=0.01)
          }
          num_per_df = nrow(chunk)
          if (!is.null(df)) {
            if (!is.null(last_epoch_st) & !is.null(last_chunk)) {
              df = clip_data(last_chunk, start_time = last_epoch_st,
                             stop_time = last_chunk[nrow(last_chunk),1])
              df = rbind(df, chunk[1:num_per_df,])
            } else {
              df = chunk[1:num_per_df,]
            }
          } else {
            df = chunk[1:num_per_df,]
          }
          result = mims_unit(df,
                             before_df = NULL,
                             after_df = NULL,
                             epoch = epoch,
                             dynamic_range = dynamic_range,
                             output_mims_per_axis = output_mims_per_axis,
                             use_gui_progress = use_gui_progress,
                             st = st,
                             use_snapshot_to_check = use_snapshot_to_check
          )
          last_epoch_st = result[nrow(result),1]
          last_chunk = chunk
          # discard the last epoch
          result = result[1:(nrow(result) - 1),]
          results[[j]] = result
          j = j + 1
        }
        else {
          if (shiny::isRunning()) {
            chunk_pb$set(value=1)
            chunk_pb$close()
          }
          close_con()
          break
        }
      }
    }
    if (shiny::isRunning()) {
      file_pb$close()
    }
    result <- do.call(rbind, results)
    return(result)
  }


.get_meta <- function(filepath, file_type, header) {
  if (file_type == 'actigraph') {
    meta = import_actigraph_meta(filepath, header = header)
    st = meta$st
    sr = meta$sr
    return(list(st = st, sr = sr))
  } else if (file_type == 'mhealth') {
    # use the first 200 rows to compute sampling rate
    funcs = import_mhealth_csv_chunked(filepath, chunk_samples = 200)
    next_chunk = funcs[[1]]
    close_con = funcs[[2]]
    df = next_chunk()
    close_con()
    sr = sampling_rate(df)
    st = df[1, 1]
    return(list(st = st, sr = sr))
  }
}