introduce detect changes api

M3nin0 · M3nin0 · commit ce544aa031ac · 2024-05-14T15:34:56.000-03:00
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -9,6 +9,10 @@ weighted_uncert_probs <- function(data_lst, unc_lst) {
     .Call(`_sits_weighted_uncert_probs`, data_lst, unc_lst)
 }
 
+dtw_distance <- function(ts1, ts2) {
+    .Call(`_sits_dtw_distance`, ts1, ts2)
+}
+
 C_kernel_median <- function(x, ncols, nrows, band, window_size) {
     .Call(`_sits_C_kernel_median`, x, ncols, nrows, band, window_size)
 }
diff --git a/R/sits_detect_change.R b/R/sits_detect_change.R
@@ -0,0 +1,60 @@
+#' @title Detect changes in time series
+#' @name sits_detect_change
+#'
+#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
+#' @author Felipe Carlos, \email{efelipecarlos@@gmail.com}
+#' @author Felipe Carvalho, \email{felipe.carvalho@@inpe.br}
+#'
+#' @description Given a set of time series or an image, this function compares
+#' each time series with a set of change/no-change patterns, and indicates
+#' places and dates where change has been detected.
+#'
+#' @param data         Set of time series.
+#' @param cd_method    Change detection method (with parameters).
+#' @param ...          Other relevant parameters.
+#' @param filter_fn    Smoothing filter function to be applied to the data.
+#' @param multicores   Number of threads to process the time series.
+#' @param progress     Show progress bar?
+#' @return             Set of time series where significant change has been
+#'                     detected.
+#' @export
+sits_detect_change <- function(data,
+                               cd_method,
+                               ...,
+                               filter_fn = NULL,
+                               multicores = 2L,
+                               progress = TRUE) {
+    UseMethod("sits_detect_change", data)
+}
+
+#' @rdname sits_detect_change
+#' @export
+sits_detect_change.sits <- function(data,
+                                    cd_method,
+                                    ...,
+                                    filter_fn = NULL,
+                                    multicores = 2L,
+                                    progress = TRUE) {
+    # set caller for error messages
+    .check_set_caller("sits_detect_change_sits")
+    # Pre-conditions
+    data <- .check_samples_ts(data)
+    .check_is_sits_model(cd_method)
+    .check_int_parameter(multicores, min = 1, max = 2048)
+    .check_progress(progress)
+    # Do detection
+    detections <- .detect_change_ts(
+        samples = data,
+        cd_method = cd_method,
+        filter_fn = filter_fn,
+        multicores = multicores,
+        progress = progress
+    )
+    return(detections)
+}
+
+#' @rdname sits_detect_change
+#' @export
+sits_detect_change.default <- function(data, cd_method, ...) {
+    stop("Input should be a sits tibble or a data cube")
+}
diff --git a/R/sits_detect_change_method.R b/R/sits_detect_change_method.R
@@ -0,0 +1,35 @@
+#' @title Create detect change method.
+#' @name sits_detect_change_method
+#'
+#' @author Gilberto Camara, \email{gilberto.camara@@inpe.br}
+#' @author Felipe Carlos, \email{efelipecarlos@@gmail.com}
+#'
+#' @description Prepare detection change method. Currently, sits supports the
+#' following methods: 'dtw' (see \code{\link[sits]{sits_dtw}})
+#'
+#' @param  samples          Time series with the training samples.
+#' @param  cd_method        Change detection method.
+#' @return                  Change detection method prepared
+#'                          to be passed to
+#'                          \code{\link[sits]{sits_detect_change}}
+#' @export
+#'
+sits_detect_change_method <- function(samples, cd_method = sits_dtw()) {
+    # set caller to show in errors
+    .check_set_caller("sits_detect_change_method")
+    # check if samples are valid
+    .check_samples_train(samples)
+    # is the train method a function?
+    .check_that(inherits(cd_method, "function"),
+                msg = .conf("messages", "sits_detect_change_method_model")
+    )
+    # are the timelines OK?
+    timeline_ok <- .timeline_check(samples)
+    .check_that(timeline_ok,
+                msg = .conf("messages", "sits_detect_change_method_timeline")
+    )
+    # compute the training method by the given data
+    result <- cd_method(samples)
+    # return a valid detect change method
+    return(result)
+}
diff --git a/man/sits_detect_change.Rd b/man/sits_detect_change.Rd
diff --git a/man/sits_detect_change_method.Rd b/man/sits_detect_change_method.Rd
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
@@ -36,6 +36,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// dtw_distance
+double dtw_distance(const NumericMatrix& ts1, const NumericMatrix& ts2);
+RcppExport SEXP _sits_dtw_distance(SEXP ts1SEXP, SEXP ts2SEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const NumericMatrix& >::type ts1(ts1SEXP);
+    Rcpp::traits::input_parameter< const NumericMatrix& >::type ts2(ts2SEXP);
+    rcpp_result_gen = Rcpp::wrap(dtw_distance(ts1, ts2));
+    return rcpp_result_gen;
+END_RCPP
+}
 // C_kernel_median
 NumericVector C_kernel_median(const NumericMatrix& x, int ncols, int nrows, int band, int window_size);
 RcppExport SEXP _sits_C_kernel_median(SEXP xSEXP, SEXP ncolsSEXP, SEXP nrowsSEXP, SEXP bandSEXP, SEXP window_sizeSEXP) {
@@ -672,6 +684,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_sits_weighted_probs", (DL_FUNC) &_sits_weighted_probs, 2},
     {"_sits_weighted_uncert_probs", (DL_FUNC) &_sits_weighted_uncert_probs, 2},
+    {"_sits_dtw_distance", (DL_FUNC) &_sits_dtw_distance, 2},
     {"_sits_C_kernel_median", (DL_FUNC) &_sits_C_kernel_median, 5},
     {"_sits_C_kernel_mean", (DL_FUNC) &_sits_C_kernel_mean, 5},
     {"_sits_C_kernel_sd", (DL_FUNC) &_sits_C_kernel_sd, 5},
diff --git a/src/detect_change_distances.cpp b/src/detect_change_distances.cpp
@@ -0,0 +1,58 @@
+#include <Rcpp.h>
+#include "./dtw.h"
+
+using namespace Rcpp;
+
+/**
+ * Convert `NumericMatrix` to 2D `std::vector`.
+ *
+ * @description
+ * This function converts a `NumericMatrix` into a 2D `std::vector`.
+ *
+ * @param mat A `NumericMatrix` with single or multi variate time-series.
+ */
+std::vector<std::vector<double>> to_cpp_vector(NumericMatrix mat) {
+    size_t rows = mat.nrow();
+    size_t cols = mat.ncol();
+
+    std::vector<std::vector<double>> result(rows, std::vector<double>(cols));
+
+    for(size_t i = 0; i < rows; ++i) {
+        for(size_t j = 0; j < cols; ++j) {
+            result[i][j] = mat(i, j);
+        }
+    }
+
+    return result;
+}
+
+/**
+ * Dynamic Time Warping (DTW) distance.
+ *
+ * @description
+ * This function calculates the Dynamic Time Warping (DTW) distance between
+ * two time-series.
+ *
+ * @param x A `double *` Time-series data.
+ * @param y A `double *` Self-Organizing Maps (SOM) codebook.
+ * @param np `int` Number of points in arrays `p1` and `p2`.
+ * @param nNA `int` Number of `NA` values in the arrays `p1` and `p2`.
+ *
+ * @reference
+ * Giorgino, T. (2009). Computing and Visualizing Dynamic Time Warping
+ * Alignments in R: The dtw Package. Journal of Statistical Software, 31(7),
+ * 1–24. https://doi.org/10.18637/jss.v031.i07
+ *
+ * @return DTW distance.
+ */
+// [[Rcpp::export]]
+double dtw_distance(
+    const NumericMatrix& ts1,
+    const NumericMatrix& ts2
+)
+{
+    std::vector<std::vector<double>> ts1_vec = to_cpp_vector(ts1);
+    std::vector<std::vector<double>> ts2_vec = to_cpp_vector(ts2);
+
+    return (distance_dtw_op(ts1_vec, ts2_vec, 2));
+}
diff --git a/src/dtw.cpp b/src/dtw.cpp
@@ -0,0 +1,100 @@
+#include <Rcpp.h>
+#include "./dtw.h"
+
+using namespace Rcpp;
+
+/**
+ * Compute the p-norm between two time-series.
+ *
+ * @description
+ * The `p-norm`, also known as the `Minkowski space`, is a generalized norm
+ * calculation that includes several types of distances based on the value
+ * of `p`.
+ *
+ * Common values of `p` include:
+ *
+ *  - `p = 1` for the Manhattan (city block) distance;
+ *  - `p = 2` for the Euclidean norm (distance).
+ *
+ * More details about p-norms can be found on Wikipedia:
+ * https://en.wikipedia.org/wiki/Norm_(mathematics)#p-norm
+ *
+ * @param a A `std::vector<double>` with time-series values.
+ * @param b A `std::vector<double>` with time-series values.
+ * @param p A `double` value of the norm to use, determining the type of
+ *          distance calculated.
+ *
+ * @note
+ * Both vectors `a` and `b` must have the same length.
+ *
+ * @note
+ * The implementation of this DTW distance calculation was adapted from the
+ * `DTW_cpp` single header library (https://github.com/cjekel/DTW_cpp).
+ *
+ * @return The `p-norm` value between vectors `a` and `b`.
+ */
+double p_norm(std::vector<double> a, std::vector<double> b, double p)
+{
+    double d = 0;
+
+    size_t index;
+    size_t a_size = a.size();
+
+    for (index = 0; index < a_size; index++)
+    {
+        d += std::pow(std::abs(a[index] - b[index]), p);
+    }
+    return std::pow(d, 1.0 / p);
+}
+
+/**
+ * Dynamic Time Warping (DTW) distance.
+ *
+ * @description
+ * This function calculates the Dynamic Time Warping (DTW) distance between
+ * two time-series.
+ *
+ * @param x A `std::vector<std::vector<double>>` with time-series values.
+ * @param y A `std::vector<std::vector<double>>` with time-series values.
+ *
+ * @reference
+ * Giorgino, T. (2009). Computing and Visualizing Dynamic Time Warping
+ * Alignments in R: The dtw Package. Journal of Statistical Software, 31(7),
+ * 1–24. https://doi.org/10.18637/jss.v031.i07
+ *
+ * @note
+ * The implementation of this DTW distance calculation was adapted from the
+ * `DTW_cpp` single header library (https://github.com/cjekel/DTW_cpp).
+ *
+ * @return DTW distance.
+ */
+double distance_dtw_op(std::vector<std::vector<double>> a,
+                       std::vector<std::vector<double>> b,
+                       double p)
+{
+    int n = a.size();
+    int o = b.size();
+
+    std::vector<std::vector<double>> d(n, std::vector<double>(o, 0.0));
+
+    d[0][0] = p_norm(a[0], b[0], p);
+
+    for (int i = 1; i < n; i++)
+    {
+        d[i][0] = d[i - 1][0] + p_norm(a[i], b[0], p);
+    }
+    for (int i = 1; i < o; i++)
+    {
+        d[0][i] = d[0][i - 1] + p_norm(a[0], b[i], p);
+    }
+    for (int i = 1; i < n; i++)
+    {
+        for (int j = 1; j < o; j++)
+        {
+            d[i][j] = p_norm(a[i], b[j], p) + std::fmin(
+                std::fmin(d[i - 1][j], d[i][j - 1]), d[i - 1][j - 1]
+            );
+        }
+    }
+    return d[n - 1][o - 1];
+}
diff --git a/src/dtw.h b/src/dtw.h
@@ -0,0 +1,5 @@
+
+#pragma once
+double distance_dtw_op(std::vector<std::vector<double>> a,
+                       std::vector<std::vector<double>> b,
+                       double p);
diff --git a/src/kohonen_distances.cpp b/src/kohonen_distances.cpp

Original file line number	Diff line number	Diff line change
`@@ -9,6 +9,10 @@ weighted_uncert_probs <- function(data_lst, unc_lst) {`
`9`	`9`	.Call(`_sits_weighted_uncert_probs`, data_lst, unc_lst)
`10`	`10`	`}`
`11`	`11`
	`12`	`+dtw_distance <- function(ts1, ts2) {`
	`13`	+ .Call(`_sits_dtw_distance`, ts1, ts2)
	`14`	`+}`
	`15`	`+`
`12`	`16`	`C_kernel_median <- function(x, ncols, nrows, band, window_size) {`
`13`	`17`	.Call(`_sits_C_kernel_median`, x, ncols, nrows, band, window_size)
`14`	`18`	`}`