snow_days_comparison.cpp

#include <unistd.h>
#include <Rcpp.h>
#include <cmath>
#include <algorithm>
#if defined(_OPENMP)
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
#endif
#include <fstream>
#include <array>
// [[Rcpp::plugins(cpp11)]]

using namespace Rcpp;


class interrupt_exception : public std::exception {
public:
  interrupt_exception(std::string message)
    : detailed_message(message)
  {};

  virtual ~interrupt_exception() throw() {};

  virtual const char* what() const throw() {
    return detailed_message.c_str();
  }

  std::string detailed_message;
};

static inline void check_interrupt_impl(void* /*dummy*/) {
  R_CheckUserInterrupt();
}

inline bool check_interrupt() {
  return (R_ToplevelExec(check_interrupt_impl, NULL) == FALSE);
}


//' xxx
//'
//' 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]]
Rcpp::List snowDaysComparisonCpp(const NumericMatrix data_1,
                                  const NumericMatrix data_2,
                                  const Rcpp::List indices,
                                  const size_t ncores,
                                  const double swe_threshold = 0.005) {

  bool interrupt = false;

  // Convert R object to proper c++ object since R objects are not thread safe
  const size_t num_pixels=data_1.nrow();
  const size_t num_years=indices.size();
  std::vector< std::vector<size_t> > indices_cpp(num_years);
  for(size_t i=0; i< num_years; ++i)
  {
    indices_cpp[i]=Rcpp::as< std::vector<size_t> >(indices[i]);
  }
  const std::vector<double> data_1_cpp = Rcpp::as< std::vector<double> >((data_1));
  const std::vector<double> data_2_cpp = Rcpp::as< std::vector<double> >((data_2));

  std::vector< std::vector< size_t > > out_cpp(num_pixels);
  for(size_t i=0; i< num_pixels; ++i){
   out_cpp[i] = std::vector< size_t >(num_years*4);
  }
  std::vector<size_t> snow_days_1(indices_cpp[0].size());
  std::vector<size_t> snow_days_2(indices_cpp[0].size());

#if defined(_OPENMP)
#pragma omp parallel num_threads(ncores) shared(interrupt, data_1_cpp, data_2_cpp, indices_cpp, out_cpp, num_pixels, num_years) private(snow_days_1, snow_days_2)
#pragma omp for schedule(static, 5)
#endif
  for(size_t i = 0; i < num_pixels; ++i){
    //Check for interupt
    if (interrupt) {
      continue;
    }
#if defined(_OPENMP)
    if (omp_get_thread_num() == 0) // only in master thread!
#endif
    if (check_interrupt()) {
      interrupt = true;
    }
    for(size_t year = 0; year < num_years; ++year) {
      snow_days_1.resize(indices_cpp[year].size());
      snow_days_2.resize(indices_cpp[year].size());
      std::fill(snow_days_1.begin(), snow_days_1.end(), 0);
      std::fill(snow_days_2.begin(), snow_days_2.end(), 0);
      size_t j = 0;
      for(const auto indice : indices_cpp[year]) {
        snow_days_1[j] = data_1_cpp[i + (indice-1)*num_pixels]>swe_threshold?1:0;
        snow_days_2[j] = data_2_cpp[i + (indice-1)*num_pixels]>swe_threshold?2:0;
        ++j;
      }
      std::transform (snow_days_1.begin(), snow_days_1.end(), snow_days_2.begin(),
                      snow_days_1.begin(), std::plus<size_t>());
      for(size_t sum = 0; sum < 4; ++sum) {
         out_cpp[i][4*year + sum] = std::count_if(snow_days_1.begin(), snow_days_1.end(),
                                                  [=](const double val){return (val == sum);});
      }
    }
  }

  if (interrupt) {
    throw interrupt_exception("[I] The computation was interrupted by user");
  }

  Rcpp::List out(num_pixels);
  Rcpp::NumericMatrix content(4,num_years);
  Rcpp::CharacterVector out_names = Rcpp::as<CharacterVector>(indices.names());
  for(size_t i = 0; i < num_pixels; ++i){
    out(i) = Rcpp::NumericMatrix(4 , num_years , out_cpp[i].begin() );
    colnames(out(i)) = out_names;
  }
  return(out);
}