store distances as double rather than float

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ripserr
 Title: Calculate Persistent Homology with Ripser-Based Engines
-Version: 1.0.0
+Version: 1.0.0.0003
 Authors@R:
     c(person(given = "Raoul R.",
              family = "Wadhwa",

NAMESPACE

@@ -2,6 +2,7 @@
 
 S3method(cubical,array)
 S3method(cubical,matrix)
+S3method(cubical,numeric)
 S3method(head,PHom)
 S3method(print,PHom)
 S3method(tail,PHom)
@@ -17,6 +18,7 @@ export(as.PHom)
 export(cubical)
 export(cubical.array)
 export(cubical.matrix)
+export(cubical.numeric)
 export(is.PHom)
 export(vietoris_rips)
 export(vietoris_rips.data.frame)

NEWS.md

@@ -1,17 +1,30 @@
 # next version
 
-## permission for deaths before births and superlevel set filtrations of rasters
+## Vietoris-Rips PH
 
-A logical argument `sublevel` has been added to `cubical` that, when `FALSE`, will pre- and post-transform raster data in order to obtain superlevel set persistent homology.
-Enabling this, an assertion that all `birth < death` has been removed from checks of persistence data.
+## float to double
+
+Ripser stores values of the `value_t` type and `ratio` as floats. This is not incompatible with R, but R users are likely to expect numeric values to be handled as doubles. Both values are now stored as doubles.
 
-## sliding window embeddings of multivariable time series (breaking change)
+### sliding window embeddings of multivariable time series (breaking change)
 
 Previously only univariable time series could be passed to `vietoris_rips()` via the sliding window embedding (used for quasi-attractor detection).
 An additional unexported embedding function has been written to handle multivariable time series. (It underperforms the original function on univariable time series, which therefore continue to rely on the original.)
 
 Furthermore, whereas `data_dim` previously defaulted to `2`, it now defaults to the number of observations per time unit of a time series as recovered by `tsp()`. (The behavior for unclassed numeric vectors remains unchanged.)
 
+## cubical PH
+
+### functionality for 1-dimensional arrays
+
+`cubical()` can now handle 1-dimensional arrays (for which no dedicated source code exists) by treating them as 2-dimensional (with an expanse of 1 in the second dimension).
+This enables the new method `cubical.numeric()` to accept vectors.
+
+### deaths before births and superlevel set filtrations
+
+A logical argument `sublevel` has been added to `cubical` that, when `FALSE`, will pre- and post-transform raster data in order to obtain superlevel set persistent homology.
+Enabling this, an assertion that all `birth < death` has been removed from checks of persistence data.
+
 # ripserr 1.0.0
 
 This major version replaces an outdated version of the Ripser C++ library with its current version.

R/cubical.R

@@ -24,6 +24,11 @@
 #' @return `PHom` object
 #' @examples 
 #' 
+#' # 1-dim example
+#' dataset <- rnorm(1500)
+#' dim(dataset) <- 1500
+#' cubical_hom1 <- cubical(dataset)
+#' 
 #' # 2-dim example
 #' dataset <- rnorm(10 ^ 2)
 #' dim(dataset) <- rep(10, 2)
@@ -74,6 +79,8 @@ cubical.array <- function(
                       method = method)
   validate_arr_cub(dataset)
 
+  # if dataset is 1-dimensional, treat it as 2-dimensional
+  if (length(dim(dataset)) == 1L) dim(dataset) <- c(dim(dataset), 1L)
 
   # transform method parameter for C++ function
   method_int <- switch(method,
@@ -146,4 +153,17 @@ cubical.matrix <- function(dataset, ...) {
 
   # return
   return(ans)
-}
+}
+
+#' @rdname cubical
+#' @export cubical.numeric
+#' @export
+cubical.numeric <- function(dataset, ...) {
+  # convert the numeric vector to a 1-dimensional array
+  dataset <- as.array(dataset, dim = 1L)
+
+  # calculate persistent homology using cubical.array
+  ans <- cubical.array(dataset, ...)
+
+  return(ans)
+}

R/utility.R

@@ -122,7 +122,7 @@ validate_arr_cub <- function(dataset) {
   error_class(dataset, "dataset", "array")
 
   # dataset should have either 2, 3, or 4 dimensions (only ones supported)
-  if (!(length(dim(dataset)) %in% c(2, 3, 4))) {
+  if (!(length(dim(dataset)) %in% seq(4))) {
     stop(paste("dataset parameter must have either 2, 3, or 4 dimensions,",
                "passed argument has", length(dim(dataset)), "dimensions"))
   }
@@ -139,7 +139,11 @@ validate_arr_cub <- function(dataset) {
   }
 
   # make sure dataset is not too large
-  if (length(dim(dataset)) == 2) {
+  if (length(dim(dataset)) == 1L) {
+    if (dim(dataset) > 2000) {
+      stop(paste("Max size for dim 1 = 2000; passed size =", dim(dataset)))
+    }
+  } else if (length(dim(dataset)) == 2) {
     if (dim(dataset)[1] > 2000 |
         dim(dataset)[2] > 1000) {
       stop(paste("Max size for dim 2 = 2000 x 1000; passed size =",

src/RcppExports.cpp

@@ -57,15 +57,15 @@ BEGIN_RCPP
 END_RCPP
 }
 // ripser_cpp_dist
-Rcpp::List ripser_cpp_dist(const Rcpp::NumericVector& dataset, int dim, double thresh, float ratio, int p);
+Rcpp::List ripser_cpp_dist(const Rcpp::NumericVector& dataset, int dim, double thresh, double ratio, int p);
 RcppExport SEXP _ripserr_ripser_cpp_dist(SEXP datasetSEXP, SEXP dimSEXP, SEXP threshSEXP, SEXP ratioSEXP, SEXP pSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type dataset(datasetSEXP);
     Rcpp::traits::input_parameter< int >::type dim(dimSEXP);
     Rcpp::traits::input_parameter< double >::type thresh(threshSEXP);
-    Rcpp::traits::input_parameter< float >::type ratio(ratioSEXP);
+    Rcpp::traits::input_parameter< double >::type ratio(ratioSEXP);
     Rcpp::traits::input_parameter< int >::type p(pSEXP);
     rcpp_result_gen = Rcpp::wrap(ripser_cpp_dist(dataset, dim, thresh, ratio, p));
     return rcpp_result_gen;

src/ripser.cpp

@@ -84,7 +84,7 @@ template <class Key> using hash = std::hash<Key>;
 
 #endif
 
-typedef float value_t;
+typedef double value_t;
 typedef int64_t index_t;
 typedef uint16_t coefficient_t;
 
@@ -413,7 +413,7 @@ template <typename DistanceMatrix> class ripser {
 	const DistanceMatrix dist;
 	const index_t n, dim_max;
 	const value_t threshold;
-	const float ratio;
+	const double ratio;
 	const coefficient_t modulus;
 	const binomial_coeff_table binomial_coeff;
 	const std::vector<coefficient_t> multiplicative_inverse;
@@ -436,7 +436,7 @@ template <typename DistanceMatrix> class ripser {
   // ripserq: Accumulate pairs in an object to be returned to the user.
   std::vector<std::vector<std::pair<value_t, value_t>>> persistence_pairs;
 
-	ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, float _ratio,
+	ripser(DistanceMatrix&& _dist, index_t _dim_max, value_t _threshold, double _ratio,
 	       coefficient_t _modulus)
 	    : dist(std::move(_dist)), n(dist.size()),
 	      dim_max(std::min(_dim_max, index_t(dist.size() - 2))), threshold(_threshold),
@@ -1257,7 +1257,7 @@ int main(int argc, char** argv) {
 
 	index_t dim_max = 1;
 	value_t threshold = std::numeric_limits<value_t>::max();
-	float ratio = 1;
+	double ratio = 1;
 	coefficient_t modulus = 2;
 
 	for (index_t i = 1; i < argc; ++i) {
@@ -1382,7 +1382,7 @@ int main(int argc, char** argv) {
 #endif
 
 // [[Rcpp::export()]]
-Rcpp::List ripser_cpp_dist(const Rcpp::NumericVector &dataset, int dim, double thresh, float ratio, int p) {
+Rcpp::List ripser_cpp_dist(const Rcpp::NumericVector &dataset, int dim, double thresh, double ratio, int p) {
   std::vector<value_t> distances(dataset.begin(), dataset.end());
 
   compressed_lower_distance_matrix dist(compressed_upper_distance_matrix(std::move(distances)));