WSL/SLF GitLab Repository

snow_days_comparison.cpp 3.82 KB
Newer Older
1
2
3
4
5
6
7
8
9
#include <unistd.h>
#include <Rcpp.h>
#include <cmath>
#include <algorithm>
#if defined(_OPENMP)
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
#endif
#include <fstream>
Adrien Michel's avatar
Adrien Michel committed
10
#include <array>
11
12
13
14
15
// [[Rcpp::plugins(cpp11)]]

using namespace Rcpp;


16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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;
};

31
32
33
34
35
36
37
38
39
static inline void check_interrupt_impl(void* /*dummy*/) {
  R_CheckUserInterrupt();
}

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


40
41
42
43
44
45
46
47
//' 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]]
48
Rcpp::List snowDaysComparisonCpp(const NumericMatrix data_1,
49
50
51
52
53
                                  const NumericMatrix data_2,
                                  const Rcpp::List indices,
                                  const size_t ncores,
                                  const double swe_threshold = 0.005) {

Adrien Michel's avatar
Adrien Michel committed
54
55
  bool interrupt = false;

56
57
58
59
60
61
62
63
64
65
  // 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));
Adrien Michel's avatar
Adrien Michel committed
66

67
68
  std::vector< std::vector< size_t > > out_cpp(num_pixels);
  for(size_t i=0; i< num_pixels; ++i){
69
   out_cpp[i] = std::vector< size_t >(num_years*4);
70
71
72
  }
  std::vector<size_t> snow_days_1(indices_cpp[0].size());
  std::vector<size_t> snow_days_2(indices_cpp[0].size());
Adrien Michel's avatar
Adrien Michel committed
73

74
75
76
77
#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
Adrien Michel's avatar
Adrien Michel committed
78
79
80
81
82
  for(size_t i = 0; i < num_pixels; ++i){
    //Check for interupt
    if (interrupt) {
      continue;
    }
83
84
85
#if defined(_OPENMP)
    if (omp_get_thread_num() == 0) // only in master thread!
#endif
Adrien Michel's avatar
Adrien Michel committed
86
87
88
89
90
91
92
93
94
95
96
97
    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;
Adrien Michel's avatar
Adrien Michel committed
98
        ++j;
Adrien Michel's avatar
Adrien Michel committed
99
100
101
      }
      std::transform (snow_days_1.begin(), snow_days_1.end(), snow_days_2.begin(),
                      snow_days_1.begin(), std::plus<size_t>());
102
      for(size_t sum = 0; sum < 4; ++sum) {
103
104
         out_cpp[i][4*year + sum] = std::count_if(snow_days_1.begin(), snow_days_1.end(),
                                                  [=](const double val){return (val == sum);});
Adrien Michel's avatar
Adrien Michel committed
105
106
107
      }
    }
  }
108
109
110
111

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

113
114
115
  Rcpp::List out(num_pixels);
  Rcpp::NumericMatrix content(4,num_years);
  Rcpp::CharacterVector out_names = Rcpp::as<CharacterVector>(indices.names());
Adrien Michel's avatar
Adrien Michel committed
116
  for(size_t i = 0; i < num_pixels; ++i){
117
118
    out(i) = Rcpp::NumericMatrix(4 , num_years , out_cpp[i].begin() );
    colnames(out(i)) = out_names;
Adrien Michel's avatar
Adrien Michel committed
119
  }
120
  return(out);
Adrien Michel's avatar
Adrien Michel committed
121
}