Changed quantile files naming

NAMESPACE

@@ -5,7 +5,6 @@ S3method(plot,vectdata)
 S3method(plot,vectdatamask)
 export(apply.model)
 export(applyQMCpp)
-export(applyQMCppDoy)
 export(assessment.plot)
 export(assessment.plot.calib)
 export(compute.data.season.mean)
@@ -13,7 +12,6 @@ export(compute.quantiles)
 export(compute.spatial.pools)
 export(compute.temporal.pools)
 export(computeQuantilesCpp)
-export(computeQuantilesCppDoy)
 export(computeSeasonnalSummaryCpp)
 export(get.filter.mask.vectdata)
 export(get.mask.na)

R/RcppExports.R

@@ -8,30 +8,8 @@
 #' everything else as FALSE.
 #' @param x An integer vector
 #' @export
-applyQMCpp <- function(swe, doyArray, sourcePath, targetPath, ncores) {
-    .Call(`_SnowQM_applyQMCpp`, swe, doyArray, sourcePath, targetPath, ncores)
-}
-
-#' applyQMCpp
-#'
-#' This function returns a logical vector identifying if
-#' there are leading NA, marking the leadings NA as TRUE and
-#' everything else as FALSE.
-#' @param x An integer vector
-#' @export
-applyQMCppDoy <- function(swe, doySelection, sourcePath, targetPath, ncores) {
-    .Call(`_SnowQM_applyQMCppDoy`, swe, doySelection, sourcePath, targetPath, ncores)
-}
-
-#' computeQuantilesCpp
-#'
-#' This function returns a logical vector identifying if
-#' there are leading NA, marking the leadings NA as TRUE and
-#' everything else as FALSE.
-#' @param x An integer vector
-#' @export
-computeQuantilesCpp <- function(swe, timePools, spatialPools, probs, outputPath, ncores) {
-    invisible(.Call(`_SnowQM_computeQuantilesCpp`, swe, timePools, spatialPools, probs, outputPath, ncores))
+applyQMCpp <- function(swe, doyArray, xCoords, yCoords, sourcePath, targetPath, ncores) {
+    .Call(`_SnowQM_applyQMCpp`, swe, doyArray, xCoords, yCoords, sourcePath, targetPath, ncores)
 }
 
 #' computeQuantilesCpp
@@ -41,8 +19,8 @@ computeQuantilesCpp <- function(swe, timePools, spatialPools, probs, outputPath,
 #' everything else as FALSE.
 #' @param x An integer vector
 #' @export
-computeQuantilesCppDoy <- function(swe, timePools, spatialPools, probs, outputPath, ncores) {
-    invisible(.Call(`_SnowQM_computeQuantilesCppDoy`, swe, timePools, spatialPools, probs, outputPath, ncores))
+computeQuantilesCpp <- function(swe, timePools, spatialPools, probs, xCoords, yCoords, outputPath, ncores) {
+    invisible(.Call(`_SnowQM_computeQuantilesCpp`, swe, timePools, spatialPools, probs, xCoords, yCoords, outputPath, ncores))
 }
 
 #' xxx

R/qm_quantile_application.R

@@ -34,6 +34,8 @@ apply.model <- function(in.data.file, out.data.file, quantile.path, ncores)
   start.time <- Sys.time()
   corrected.swe <- applyQMCpp(data$data,
                               doyArray = doys,
+                              xCoords=round(data$x.coords[data$x.indices]),
+                              yCoords=round(data$y.coords[data$y.indices]),
                               sourcePath = paste0(quantile.path,"/training_quantiles"),
                               targetPath = paste0(quantile.path,"/target_quantiles"),
                               ncores = ncores)

R/qm_quantile_computation.R

@@ -252,9 +252,9 @@ compute.quantiles <- function(data, temporal.pools, spatial.pools, output.path,
   computeQuantilesCpp(swe=data$data,
                       timePools=temporal.pools,
                       spatialPools=spatial.pools,
-                      xCoords=round(data$x.coords[data$x.indices]),
-                      yCoords=round(data$x.yoords[data$xy.indices]),
                       probs=probs,
+                      xCoords=round(data$x.coords[data$x.indices]),
+                      yCoords=round(data$y.coords[data$y.indices]),
                       outputPath=output.path, ncores=ncores)
 
 }

man/applyQMCpp.Rd

@@ -4,7 +4,7 @@
 \alias{applyQMCpp}
 \title{applyQMCpp}
 \usage{
-applyQMCpp(swe, doyArray, sourcePath, targetPath, ncores)
+applyQMCpp(swe, doyArray, xCoords, yCoords, sourcePath, targetPath, ncores)
 }
 \arguments{
 \item{x}{An integer vector}

man/applyQMCppDoy.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
-\name{applyQMCppDoy}
-\alias{applyQMCppDoy}
-\title{applyQMCpp}
-\usage{
-applyQMCppDoy(swe, doySelection, sourcePath, targetPath, ncores)
-}
-\arguments{
-\item{x}{An integer vector}
-}
-\description{
-This function returns a logical vector identifying if
-there are leading NA, marking the leadings NA as TRUE and
-everything else as FALSE.
-}

man/computeQuantilesCpp.Rd

@@ -4,7 +4,16 @@
 \alias{computeQuantilesCpp}
 \title{computeQuantilesCpp}
 \usage{
-computeQuantilesCpp(swe, timePools, spatialPools, probs, outputPath, ncores)
+computeQuantilesCpp(
+  swe,
+  timePools,
+  spatialPools,
+  probs,
+  xCoords,
+  yCoords,
+  outputPath,
+  ncores
+)
 }
 \arguments{
 \item{x}{An integer vector}

man/computeQuantilesCppDoy.Rd

-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RcppExports.R
-\name{computeQuantilesCppDoy}
-\alias{computeQuantilesCppDoy}
-\title{computeQuantilesCpp}
-\usage{
-computeQuantilesCppDoy(swe, timePools, spatialPools, probs, outputPath, ncores)
-}
-\arguments{
-\item{x}{An integer vector}
-}
-\description{
-This function returns a logical vector identifying if
-there are leading NA, marking the leadings NA as TRUE and
-everything else as FALSE.
-}

man/read.mask.asc.Rd

@@ -4,7 +4,7 @@
 \alias{read.mask.asc}
 \title{Title}
 \usage{
-read.mask.asc(filename, val.to.remove = NULL, x.shift = NULL, y.shift = NULL)
+read.mask.asc(filename, val.to.remove = NULL, x.shift = 0, y.shift = 0)
 }
 \arguments{
 \item{filename}{A}

man/read.mask.tif.Rd

@@ -4,7 +4,7 @@
 \alias{read.mask.tif}
 \title{Title}
 \usage{
-read.mask.tif(filename, val.to.remove = NULL, x.shift = NULL, y.shift = NULL)
+read.mask.tif(filename, val.to.remove = NULL, x.shift = 0, y.shift = 0)
 }
 \arguments{
 \item{filename}{A}

man/read.swe.Rd

@@ -16,6 +16,7 @@ read.swe(
   dem.path,
   slope.path,
   azi.path,
+  curvature.path,
   x.shift = 0,
   y.shift = 0,
   silent = FALSE

src/RcppExports.cpp

@@ -11,62 +11,36 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
 // applyQMCpp
-NumericVector applyQMCpp(const Rcpp::NumericMatrix swe, const Rcpp::NumericVector doyArray, const std::string sourcePath, const std::string targetPath, const size_t ncores);
-RcppExport SEXP _SnowQM_applyQMCpp(SEXP sweSEXP, SEXP doyArraySEXP, SEXP sourcePathSEXP, SEXP targetPathSEXP, SEXP ncoresSEXP) {
+NumericVector applyQMCpp(const Rcpp::NumericMatrix swe, const Rcpp::NumericVector doyArray, const NumericVector xCoords, const NumericVector yCoords, const std::string sourcePath, const std::string targetPath, const size_t ncores);
+RcppExport SEXP _SnowQM_applyQMCpp(SEXP sweSEXP, SEXP doyArraySEXP, SEXP xCoordsSEXP, SEXP yCoordsSEXP, SEXP sourcePathSEXP, SEXP targetPathSEXP, SEXP ncoresSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< const Rcpp::NumericMatrix >::type swe(sweSEXP);
     Rcpp::traits::input_parameter< const Rcpp::NumericVector >::type doyArray(doyArraySEXP);
+    Rcpp::traits::input_parameter< const NumericVector >::type xCoords(xCoordsSEXP);
+    Rcpp::traits::input_parameter< const NumericVector >::type yCoords(yCoordsSEXP);
     Rcpp::traits::input_parameter< const std::string >::type sourcePath(sourcePathSEXP);
     Rcpp::traits::input_parameter< const std::string >::type targetPath(targetPathSEXP);
     Rcpp::traits::input_parameter< const size_t >::type ncores(ncoresSEXP);
-    rcpp_result_gen = Rcpp::wrap(applyQMCpp(swe, doyArray, sourcePath, targetPath, ncores));
-    return rcpp_result_gen;
-END_RCPP
-}
-// applyQMCppDoy
-NumericVector applyQMCppDoy(const Rcpp::NumericMatrix swe, const Rcpp::List doySelection, const std::string sourcePath, const std::string targetPath, const size_t ncores);
-RcppExport SEXP _SnowQM_applyQMCppDoy(SEXP sweSEXP, SEXP doySelectionSEXP, SEXP sourcePathSEXP, SEXP targetPathSEXP, SEXP ncoresSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const Rcpp::NumericMatrix >::type swe(sweSEXP);
-    Rcpp::traits::input_parameter< const Rcpp::List >::type doySelection(doySelectionSEXP);
-    Rcpp::traits::input_parameter< const std::string >::type sourcePath(sourcePathSEXP);
-    Rcpp::traits::input_parameter< const std::string >::type targetPath(targetPathSEXP);
-    Rcpp::traits::input_parameter< const size_t >::type ncores(ncoresSEXP);
-    rcpp_result_gen = Rcpp::wrap(applyQMCppDoy(swe, doySelection, sourcePath, targetPath, ncores));
+    rcpp_result_gen = Rcpp::wrap(applyQMCpp(swe, doyArray, xCoords, yCoords, sourcePath, targetPath, ncores));
     return rcpp_result_gen;
 END_RCPP
 }
 // computeQuantilesCpp
-void computeQuantilesCpp(const NumericMatrix swe, const Rcpp::List timePools, const Rcpp::List spatialPools, const NumericVector probs, const std::string outputPath, const size_t ncores);
-RcppExport SEXP _SnowQM_computeQuantilesCpp(SEXP sweSEXP, SEXP timePoolsSEXP, SEXP spatialPoolsSEXP, SEXP probsSEXP, SEXP outputPathSEXP, SEXP ncoresSEXP) {
+void computeQuantilesCpp(const Rcpp::NumericMatrix swe, const Rcpp::List timePools, const Rcpp::List spatialPools, const Rcpp::NumericVector probs, const Rcpp::NumericVector xCoords, const Rcpp::NumericVector yCoords, const std::string outputPath, const size_t ncores);
+RcppExport SEXP _SnowQM_computeQuantilesCpp(SEXP sweSEXP, SEXP timePoolsSEXP, SEXP spatialPoolsSEXP, SEXP probsSEXP, SEXP xCoordsSEXP, SEXP yCoordsSEXP, SEXP outputPathSEXP, SEXP ncoresSEXP) {
 BEGIN_RCPP
     Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const NumericMatrix >::type swe(sweSEXP);
-    Rcpp::traits::input_parameter< const Rcpp::List >::type timePools(timePoolsSEXP);
-    Rcpp::traits::input_parameter< const Rcpp::List >::type spatialPools(spatialPoolsSEXP);
-    Rcpp::traits::input_parameter< const NumericVector >::type probs(probsSEXP);
-    Rcpp::traits::input_parameter< const std::string >::type outputPath(outputPathSEXP);
-    Rcpp::traits::input_parameter< const size_t >::type ncores(ncoresSEXP);
-    computeQuantilesCpp(swe, timePools, spatialPools, probs, outputPath, ncores);
-    return R_NilValue;
-END_RCPP
-}
-// computeQuantilesCppDoy
-void computeQuantilesCppDoy(const NumericMatrix swe, const Rcpp::List timePools, const Rcpp::List spatialPools, const NumericVector probs, const std::string outputPath, const size_t ncores);
-RcppExport SEXP _SnowQM_computeQuantilesCppDoy(SEXP sweSEXP, SEXP timePoolsSEXP, SEXP spatialPoolsSEXP, SEXP probsSEXP, SEXP outputPathSEXP, SEXP ncoresSEXP) {
-BEGIN_RCPP
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const NumericMatrix >::type swe(sweSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericMatrix >::type swe(sweSEXP);
     Rcpp::traits::input_parameter< const Rcpp::List >::type timePools(timePoolsSEXP);
     Rcpp::traits::input_parameter< const Rcpp::List >::type spatialPools(spatialPoolsSEXP);
-    Rcpp::traits::input_parameter< const NumericVector >::type probs(probsSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericVector >::type probs(probsSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericVector >::type xCoords(xCoordsSEXP);
+    Rcpp::traits::input_parameter< const Rcpp::NumericVector >::type yCoords(yCoordsSEXP);
     Rcpp::traits::input_parameter< const std::string >::type outputPath(outputPathSEXP);
     Rcpp::traits::input_parameter< const size_t >::type ncores(ncoresSEXP);
-    computeQuantilesCppDoy(swe, timePools, spatialPools, probs, outputPath, ncores);
+    computeQuantilesCpp(swe, timePools, spatialPools, probs, xCoords, yCoords, outputPath, ncores);
     return R_NilValue;
 END_RCPP
 }
@@ -116,10 +90,8 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_SnowQM_applyQMCpp", (DL_FUNC) &_SnowQM_applyQMCpp, 5},
-    {"_SnowQM_applyQMCppDoy", (DL_FUNC) &_SnowQM_applyQMCppDoy, 5},
-    {"_SnowQM_computeQuantilesCpp", (DL_FUNC) &_SnowQM_computeQuantilesCpp, 6},
-    {"_SnowQM_computeQuantilesCppDoy", (DL_FUNC) &_SnowQM_computeQuantilesCppDoy, 6},
+    {"_SnowQM_applyQMCpp", (DL_FUNC) &_SnowQM_applyQMCpp, 7},
+    {"_SnowQM_computeQuantilesCpp", (DL_FUNC) &_SnowQM_computeQuantilesCpp, 8},
     {"_SnowQM_computeSeasonnalSummaryCpp", (DL_FUNC) &_SnowQM_computeSeasonnalSummaryCpp, 4},
     {"_SnowQM_snowDaysComparisonCpp", (DL_FUNC) &_SnowQM_snowDaysComparisonCpp, 5},
     {"_SnowQM_snowSeasonLengthCpp", (DL_FUNC) &_SnowQM_snowSeasonLengthCpp, 5},

src/quantile_application.cpp

@@ -90,6 +90,8 @@ bool getFileContent(const std::string fileName,
 // [[Rcpp::export]]
 NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
                          const Rcpp::NumericVector doyArray,
+                         const NumericVector xCoords,
+                         const NumericVector yCoords,
                          const std::string sourcePath,
                          const std::string targetPath,
                          const size_t ncores)
@@ -97,6 +99,8 @@ NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
   // Convert R object to proper c++ object since R objects are not thread safe
   const size_t num_doys=366;
   std::vector<size_t> doyArray_cpp = Rcpp::as< std::vector<size_t> >(doyArray);
+  const std::vector<size_t> xCoords_cpp =  Rcpp::as< std::vector<size_t> >(xCoords);
+  const std::vector<size_t> yCoords_cpp =  Rcpp::as< std::vector<size_t> >(yCoords);
 
   const size_t num_pixels=swe.nrow();
   std::vector< size_t > pixels(num_pixels);
@@ -127,13 +131,15 @@ NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
     }
     // Read QM files
     std::vector< std::vector <double>> sourcePDF;
-    if(!getFileContent(sourcePath+"/quantiles_"+std::to_string(p)+".txt",
+    if(!getFileContent(sourcePath+"/quantiles_"+std::to_string(xCoords_cpp[p])+
+                       "_"+std::to_string(yCoords_cpp[p])+".txt",
                        sourcePDF, num_doys)){
       file_io_error = true;
       continue;
     }
     std::vector< std::vector <double>> targetPDF;
-    if(!getFileContent(targetPath+"/quantiles_"+std::to_string(p)+".txt",
+    if(!getFileContent(targetPath+"/quantiles_"+std::to_string(xCoords_cpp[p])+
+                       "_"+std::to_string(yCoords_cpp[p])+".txt",
                        targetPDF, num_doys)){
       file_io_error = true;
       continue;
@@ -176,94 +182,94 @@ NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
   return wrap(swe_cpp);
 }
 
-//' applyQMCpp
-//'
-//' This function returns a logical vector identifying if
-//' there are leading NA, marking the leadings NA as TRUE and
-//' everything else as FALSE.
-//' @param x An integer vector
-//' @export
-// [[Rcpp::export]]
-NumericVector applyQMCppDoy(const Rcpp::NumericMatrix swe,
-                         const Rcpp::List doySelection,
-                         const std::string sourcePath,
-                         const std::string targetPath,
-                         const size_t ncores)
-{
-  // Convert R object to proper c++ object since R objects are not thread safe
-  const size_t num_doys=doySelection.size();
-
-  const Rcpp::CharacterVector doyRealNumber = doySelection.names();
-  std::vector< std::string > doyRealNumber_cpp(num_doys);
-  std::vector< std::vector<size_t> > doySelection_cpp(num_doys);
-  for(size_t i=0; i< num_doys;i++)
-  {
-    doyRealNumber_cpp[i] = doyRealNumber[i];
-    doySelection_cpp[i] = Rcpp::as< std::vector<size_t> >(doySelection[i]);
-  }
-  const size_t num_pixels=swe.nrow();
-  std::vector< size_t > pixels(num_pixels);
-  for(size_t i=0; i < num_pixels; i++) {
-    pixels[i]=i;
-  }
-
-  std::vector<double> swe_cpp = Rcpp::as< std::vector<double> >(swe);
-
-  bool interrupt = false;
-  bool file_io_error = false;
-#if defined(_OPENMP)
-#pragma omp parallel num_threads(ncores) shared(interrupt, pixels, num_doys, num_pixels, swe_cpp, doyRealNumber_cpp, doySelection_cpp, sourcePath, targetPath)
-#pragma omp for schedule(static, 5)
-#endif
-  for(size_t doy = 0; doy < num_doys; doy++) {
-    //Check for interupt
-    if (interrupt | file_io_error) {
-      continue;
-    }
-#if defined(_OPENMP)
-    if (omp_get_thread_num() == 0) // only in master thread!
-#endif
-    if (check_interrupt()) {
-      interrupt = true;
-      continue;
-    }
-
-    const std::string doyNumber = doyRealNumber_cpp[doy];
-    // Read QM files
-    std::vector< std::vector <double>> sourcePDF;
-    if(!getFileContent(sourcePath+"/quantiles_"+doyNumber+".txt",
-                       sourcePDF, num_pixels)){
-      file_io_error = true;
-      continue;
-    }
-    std::vector< std::vector <double>> targetPDF;
-    if(!getFileContent(targetPath+"/quantiles_"+doyNumber+".txt",
-                       targetPDF, num_pixels)){
-      file_io_error = true;
-      continue;
-    }
-    for(const auto p:pixels) {
-      if(sourcePDF[p].size() != targetPDF[p].size()) {
-        file_io_error = true;
-        continue;
-      }
-      for(const auto d:doySelection_cpp[doy]){
-        double& target = swe_cpp[p+(d-1)*num_pixels];
-        auto lower = std::upper_bound(sourcePDF[p].begin(), sourcePDF[p].end(), target);
-        size_t quantile = std::distance(sourcePDF[p].begin(),lower);
-        //Apply the -1 correction and lower bound
-        // Here quantiles are 0...100, in R 1...101
-        quantile = (quantile <= 1) ? 1 : quantile - 1;
-        quantile = (quantile >= (sourcePDF[p].size() - 1)) ? sourcePDF[p].size() - 2 : quantile;
-        target += targetPDF[p][quantile] - sourcePDF[p][quantile];
-      }
-    }
-  }
-  if (interrupt) {
-    throw interrupt_exception("[E] The computation was interrupted by user");
-  }
-  if(file_io_error){
-    throw interrupt_exception("[E] Problem reading QM files, either files do not exist, they do not match the right dimensions, or they are ill formated");
-  }
-  return wrap(swe_cpp);
-}
+// //' applyQMCpp
+// //'
+// //' This function returns a logical vector identifying if
+// //' there are leading NA, marking the leadings NA as TRUE and
+// //' everything else as FALSE.
+// //' @param x An integer vector
+// //' @export
+// // [[Rcpp::export]]
+// NumericVector applyQMCppDoy(const Rcpp::NumericMatrix swe,
+//                          const Rcpp::List doySelection,
+//                          const std::string sourcePath,
+//                          const std::string targetPath,
+//                          const size_t ncores)
+// {
+//   // Convert R object to proper c++ object since R objects are not thread safe
+//   const size_t num_doys=doySelection.size();
+//
+//   const Rcpp::CharacterVector doyRealNumber = doySelection.names();
+//   std::vector< std::string > doyRealNumber_cpp(num_doys);
+//   std::vector< std::vector<size_t> > doySelection_cpp(num_doys);
+//   for(size_t i=0; i< num_doys;i++)
+//   {
+//     doyRealNumber_cpp[i] = doyRealNumber[i];
+//     doySelection_cpp[i] = Rcpp::as< std::vector<size_t> >(doySelection[i]);
+//   }
+//   const size_t num_pixels=swe.nrow();
+//   std::vector< size_t > pixels(num_pixels);
+//   for(size_t i=0; i < num_pixels; i++) {
+//     pixels[i]=i;
+//   }
+//
+//   std::vector<double> swe_cpp = Rcpp::as< std::vector<double> >(swe);
+//
+//   bool interrupt = false;
+//   bool file_io_error = false;
+// #if defined(_OPENMP)
+// #pragma omp parallel num_threads(ncores) shared(interrupt, pixels, num_doys, num_pixels, swe_cpp, doyRealNumber_cpp, doySelection_cpp, sourcePath, targetPath)
+// #pragma omp for schedule(static, 5)
+// #endif
+//   for(size_t doy = 0; doy < num_doys; doy++) {
+//     //Check for interupt
+//     if (interrupt | file_io_error) {
+//       continue;
+//     }
+// #if defined(_OPENMP)
+//     if (omp_get_thread_num() == 0) // only in master thread!
+// #endif
+//     if (check_interrupt()) {
+//       interrupt = true;
+//       continue;
+//     }
+//
+//     const std::string doyNumber = doyRealNumber_cpp[doy];
+//     // Read QM files
+//     std::vector< std::vector <double>> sourcePDF;
+//     if(!getFileContent(sourcePath+"/quantiles_"+doyNumber+".txt",
+//                        sourcePDF, num_pixels)){
+//       file_io_error = true;
+//       continue;
+//     }
+//     std::vector< std::vector <double>> targetPDF;
+//     if(!getFileContent(targetPath+"/quantiles_"+doyNumber+".txt",
+//                        targetPDF, num_pixels)){
+//       file_io_error = true;
+//       continue;
+//     }
+//     for(const auto p:pixels) {
+//       if(sourcePDF[p].size() != targetPDF[p].size()) {
+//         file_io_error = true;
+//         continue;
+//       }
+//       for(const auto d:doySelection_cpp[doy]){
+//         double& target = swe_cpp[p+(d-1)*num_pixels];
+//         auto lower = std::upper_bound(sourcePDF[p].begin(), sourcePDF[p].end(), target);
+//         size_t quantile = std::distance(sourcePDF[p].begin(),lower);
+//         //Apply the -1 correction and lower bound
+//         // Here quantiles are 0...100, in R 1...101
+//         quantile = (quantile <= 1) ? 1 : quantile - 1;
+//         quantile = (quantile >= (sourcePDF[p].size() - 1)) ? sourcePDF[p].size() - 2 : quantile;
+//         target += targetPDF[p][quantile] - sourcePDF[p][quantile];
+//       }
+//     }
+//   }
+//   if (interrupt) {
+//     throw interrupt_exception("[E] The computation was interrupted by user");
+//   }
+//   if(file_io_error){
+//     throw interrupt_exception("[E] Problem reading QM files, either files do not exist, they do not match the right dimensions, or they are ill formated");
+//   }
+//   return wrap(swe_cpp);
+// }

src/quantile_computation.cpp

@@ -7,8 +7,6 @@
 #endif
 #include <fstream>
 
-using namespace Rcpp;
-
 class interrupt_exception : public std::exception {
 public:
   interrupt_exception(std::string message)
@@ -90,15 +88,19 @@ inline bool check_interrupt() {
 //' @param x An integer vector
 //' @export
 // [[Rcpp::export]]
-void computeQuantilesCpp(const NumericMatrix swe,
+void computeQuantilesCpp(const Rcpp::NumericMatrix swe,
                              const Rcpp::List timePools,
                              const Rcpp::List spatialPools,
-                             const NumericVector probs,
+                             const Rcpp::NumericVector probs,
+                             const Rcpp::NumericVector xCoords,
+                             const Rcpp::NumericVector yCoords,
                              const std::string outputPath,
                              const size_t ncores)
 {
   // Convert R object to proper c++ object since R objects are not thread safe
-  const std::vector<double> probs_cpp = as<std::vector<double> >(probs);
+  const std::vector<size_t> xCoords_cpp =  Rcpp::as< std::vector<size_t> >(xCoords);
+  const std::vector<size_t> yCoords_cpp =  Rcpp::as< std::vector<size_t> >(yCoords);
+  const std::vector<double> probs_cpp = Rcpp::as<std::vector<double> >(probs);
   const size_t num_pixels=swe.nrow();
   std::vector< std::vector<size_t> > spatialPools_cpp(num_pixels);
   for(size_t i=0; i< num_pixels;i++)
@@ -119,7 +121,7 @@ void computeQuantilesCpp(const NumericMatrix swe,
 #pragma omp for schedule(dynamic, 10)
 #endif
   for(size_t pixel=0; pixel<num_pixels; pixel++) {
-    
+
     //Check for interupt
     if (interrupt) {
       continue;
@@ -130,7 +132,7 @@ void computeQuantilesCpp(const NumericMatrix swe,
     if (check_interrupt()) {
       interrupt = true;
     }
-    
+
     const std::vector<size_t>& pixels_indices = spatialPools_cpp[pixel];
     const size_t num_pooled_pixels = pixels_indices.size();
 
@@ -140,7 +142,7 @@ void computeQuantilesCpp(const NumericMatrix swe,
 
       const std::vector<size_t>& time_indices = timePools_cpp[doy];
       const size_t num_pooled_timestep = time_indices.size();
-      
+
       std::vector<double> quantiles_data(num_pooled_pixels*num_pooled_timestep);
       for (size_t i=0; i < num_pooled_pixels; i++) {
         for (size_t j=0; j < num_pooled_timestep; j++) {
@@ -152,7 +154,8 @@ void computeQuantilesCpp(const NumericMatrix swe,
       out_cpp[doy] = quantile_cpp(quantiles_data,probs_cpp);
     }
 
-    std::ofstream outFile(outputPath+"/quantiles_"+std::to_string(pixel)+".txt");
+    std::ofstream outFile(outputPath+"/quantiles_"+std::to_string(xCoords_cpp[pixel])+
+                          "_"+std::to_string(yCoords_cpp[pixel])+".txt");
     for (const auto &line : out_cpp) {
       for (const auto &e : line) {
         outFile << e << " ";
@@ -165,84 +168,84 @@ void computeQuantilesCpp(const NumericMatrix swe,
   }
 }
 
-//' computeQuantilesCpp
-//'
-//' This function returns a logical vector identifying if
-//' there are leading NA, marking the leadings NA as TRUE and
-//' everything else as FALSE.
-//' @param x An integer vector
-//' @export
-// [[Rcpp::export]]
-void computeQuantilesCppDoy(const NumericMatrix swe,
-                         const Rcpp::List timePools,
-                         const Rcpp::List spatialPools,
-                         const NumericVector probs,
-                         const std::string outputPath,
-                         const size_t ncores)
-{
-  // Convert R object to proper c++ object since R objects are not thread safe
-  const std::vector<double> probs_cpp = as<std::vector<double> >(probs);
-  const size_t num_pixels=swe.nrow();
-  std::vector< std::vector<size_t> > spatialPools_cpp(num_pixels);
-  for(size_t i=0; i< num_pixels;i++)
-  {
-    spatialPools_cpp[i]=Rcpp::as< std::vector<size_t> >(spatialPools[i]);
-  }
-  std::vector< std::vector<size_t> > timePools_cpp(366);
-  for(size_t i=0; i< 366;i++)
-  {
-    timePools_cpp[i]=Rcpp::as< std::vector<size_t> >(timePools[i]);
-  }
-  const std::vector<double> swe_cpp = Rcpp::as< std::vector<double> >(swe);
-  
-  bool interrupt = false;
-  
-#if defined(_OPENMP)
-#pragma omp parallel num_threads(ncores) shared(interrupt, probs_cpp, num_pixels, spatialPools_cpp, timePools_cpp, swe_cpp)
-#pragma omp for schedule(static, 5)
-#endif
-  for(size_t doy = 1; doy <= 366; doy++) {
-    //Check for interupt
-    if (interrupt) {
-      continue;
-    }
-#if defined(_OPENMP)
-    if (omp_get_thread_num() == 0) // only in master thread!