WSL/SLF GitLab Repository

quantile_application.cpp 4.93 KB
Newer Older
Adrien Michel's avatar
Adrien Michel committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#include <unistd.h>
#include <Rcpp.h>
#include <cmath>
#include <algorithm>
#if defined(_OPENMP)
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
#endif
#include <fstream>

// [[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);
}

bool getFileContent(const std::string fileName,
                    std::vector< std::vector <double>>&  matrix,
                    size_t nrows)
{
  // Open the File
  std::ifstream file(fileName.c_str());
  if(!file)
  {
    return false;
  }

  // Check number of lines
  const size_t row_counts =  std::count(std::istreambuf_iterator<char>(file),
                             std::istreambuf_iterator<char>(), '\n');
  if(row_counts != nrows) {
    return false;
  }
  // Returns to the beginning of the file
  file.clear();
  file.seekg(0);

  // Check number of columns
  std::string line;
  std::getline(file, line);
  std::stringstream s;
  s << line;
  size_t ncols = 0;
  double summy_value;
  while(s >> summy_value) ++ncols;

  // Returns to the beginning of the file
  file.clear();
  file.seekg(0);

  matrix.resize(nrows);
  for(auto& m : matrix)
    m.resize(ncols);
  for(auto& outer : matrix)
    for(auto& inner : outer)
      file >> inner;

  file.close();
  return true;
}


//' 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 applyQMCpp(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;
139
    if(!getFileContent(sourcePath+"/quantiles_"+doyNumber+".txt",
Adrien Michel's avatar
Adrien Michel committed
140
141
142
143
144
                       sourcePDF, num_pixels)){
      file_io_error = true;
      continue;
    }
    std::vector< std::vector <double>> targetPDF;
145
    if(!getFileContent(targetPath+"/quantiles_"+doyNumber+".txt",
Adrien Michel's avatar
Adrien Michel committed
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
                       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);
}