WSL/SLF GitLab Repository

Commit 6f42e33d authored by Adrien Michel's avatar Adrien Michel
Browse files

changed loop order in quantile computation and application, silenced warnings

parent 9a6f0303
......@@ -5,6 +5,7 @@ S3method(plot,vectdata)
S3method(plot,vectdatamask)
export(apply.model)
export(applyQMCpp)
export(applyQMCppDoy)
export(assessment.plot)
export(assessment.plot.calib)
export(compute.data.season.diff.calib)
......@@ -12,6 +13,7 @@ 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)
......
......@@ -8,8 +8,19 @@
#' everything else as FALSE.
#' @param x An integer vector
#' @export
applyQMCpp <- function(swe, doySelection, sourcePath, targetPath, ncores) {
.Call(`_SnowQM_applyQMCpp`, swe, doySelection, sourcePath, targetPath, ncores)
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
......@@ -23,6 +34,17 @@ computeQuantilesCpp <- function(swe, timePools, spatialPools, probs, outputPath,
invisible(.Call(`_SnowQM_computeQuantilesCpp`, swe, timePools, spatialPools, probs, outputPath, 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
computeQuantilesCppDoy <- function(swe, timePools, spatialPools, probs, outputPath, ncores) {
invisible(.Call(`_SnowQM_computeQuantilesCppDoy`, swe, timePools, spatialPools, probs, outputPath, ncores))
}
#' xxx
#'
#' This function returns a logical vector identifying if
......
......@@ -91,7 +91,7 @@ apply.model <- function(in.data.file, out.data.file, quantile.path, type, ncores
start.time <- Sys.time()
if(type == "cpp") {
corrected.swe <- applyQMCpp(data$data,
doySelection = doy.selection,
doyArray = doys,
sourcePath = paste0(quantile.path,"/training_quantiles"),
targetPath = paste0(quantile.path,"/target_quantiles"),
ncores = ncores)
......
......@@ -4,7 +4,7 @@
\alias{applyQMCpp}
\title{applyQMCpp}
\usage{
applyQMCpp(swe, doySelection, sourcePath, targetPath, ncores)
applyQMCpp(swe, doyArray, sourcePath, targetPath, ncores)
}
\arguments{
\item{x}{An integer vector}
......
PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS)
PKG_LIBS = $(SHLIB_OPENMP_CFLAGS)
PKG_CXXFLAGS= $(SHLIB_OPENMP_CFLAGS)
PKG_CXXFLAGS= $(SHLIB_OPENMP_CFLAGS) -Wpedantic -O3
CXX_STD = CXX11
\ No newline at end of file
......@@ -11,8 +11,23 @@ Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
#endif
// applyQMCpp
NumericVector applyQMCpp(const Rcpp::NumericMatrix swe, const Rcpp::List doySelection, const std::string sourcePath, const std::string targetPath, const size_t ncores);
RcppExport SEXP _SnowQM_applyQMCpp(SEXP sweSEXP, SEXP doySelectionSEXP, SEXP sourcePathSEXP, SEXP targetPathSEXP, SEXP ncoresSEXP) {
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) {
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 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;
......@@ -21,7 +36,7 @@ BEGIN_RCPP
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, doySelection, sourcePath, targetPath, ncores));
rcpp_result_gen = Rcpp::wrap(applyQMCppDoy(swe, doySelection, sourcePath, targetPath, ncores));
return rcpp_result_gen;
END_RCPP
}
......@@ -40,6 +55,21 @@ BEGIN_RCPP
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::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);
computeQuantilesCppDoy(swe, timePools, spatialPools, probs, outputPath, ncores);
return R_NilValue;
END_RCPP
}
// computeSeasonnalSummaryCpp
Rcpp::NumericMatrix computeSeasonnalSummaryCpp(const Rcpp::NumericMatrix swe, const std::string fun, const Rcpp::List indices, const double swe_threshold);
RcppExport SEXP _SnowQM_computeSeasonnalSummaryCpp(SEXP sweSEXP, SEXP funSEXP, SEXP indicesSEXP, SEXP swe_thresholdSEXP) {
......@@ -87,7 +117,9 @@ 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_computeSeasonnalSummaryCpp", (DL_FUNC) &_SnowQM_computeSeasonnalSummaryCpp, 4},
{"_SnowQM_snowDaysComparisonCpp", (DL_FUNC) &_SnowQM_snowDaysComparisonCpp, 5},
{"_SnowQM_snowSeasonLengthCpp", (DL_FUNC) &_SnowQM_snowSeasonLengthCpp, 5},
......
......@@ -68,7 +68,6 @@ bool getFileContent(const std::string fileName,
// Returns to the beginning of the file
file.clear();
file.seekg(0);
matrix.resize(nrows);
for(auto& m : matrix)
m.resize(ncols);
......@@ -90,6 +89,92 @@ bool getFileContent(const std::string fileName,
//' @export
// [[Rcpp::export]]
NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
const Rcpp::NumericVector doyArray,
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=Rcpp::max(doyArray);
std::vector<size_t> doyArray_cpp = Rcpp::as< std::vector<size_t> >(doyArray);
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, doyArray_cpp, sourcePath, targetPath)
#pragma omp for schedule(static, 5)
#endif
for(const auto p:pixels) {
//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;
}
// Read QM files
std::vector< std::vector <double>> sourcePDF;
if(!getFileContent(sourcePath+"/quantiles_"+std::to_string(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",
targetPDF, num_doys)){
file_io_error = true;
continue;
}
for(size_t d = 0; d < num_doys; ++d){
if(sourcePDF[d].size() != targetPDF[d].size()) {
file_io_error = true;
continue;
}
}
for(size_t d = 0; d < doyArray_cpp.size(); ++d) {
const size_t doy(doyArray_cpp[d]-1);
double& target = swe_cpp[p+(d)*num_pixels];
auto lower = std::upper_bound(sourcePDF[doy].begin(), sourcePDF[doy].end(), target);
size_t quantile = std::distance(sourcePDF[doy].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[doy].size() - 1)) ? sourcePDF[doy].size() - 2 : quantile;
target += targetPDF[doy][quantile] - sourcePDF[doy][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,
......@@ -97,7 +182,7 @@ 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=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);
......@@ -111,14 +196,14 @@ NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
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)
#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
......@@ -132,7 +217,7 @@ NumericVector applyQMCpp(const Rcpp::NumericMatrix swe,
interrupt = true;
continue;
}
const std::string doyNumber = doyRealNumber_cpp[doy];
// Read QM files
std::vector< std::vector <double>> sourcePDF;
......
......@@ -3,13 +3,10 @@
#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 {
......@@ -117,6 +114,89 @@ void computeQuantilesCpp(const NumericMatrix 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(dynamic, 10)
#endif
for(size_t pixel=0; pixel<num_pixels; pixel++) {
//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;
}
const std::vector<size_t>& pixels_indices = spatialPools_cpp[pixel];
const size_t num_pooled_pixels = pixels_indices.size();
std::vector<std::vector<double>> out_cpp(366);
for(size_t doy = 0; doy < 366; doy++) {
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++) {
quantiles_data[i+j*num_pooled_pixels] = swe_cpp[(pixels_indices[i]-1) +
(time_indices[j]-1)*num_pixels];
}
}
remove_nan(quantiles_data);
out_cpp[doy] = quantile_cpp(quantiles_data,probs_cpp);
}
std::ofstream outFile(outputPath+"/quantiles_"+std::to_string(pixel)+".txt");
for (const auto &line : out_cpp) {
for (const auto &e : line) {
outFile << e << " ";
}
outFile << "\n";
}
}
if (interrupt) {
throw interrupt_exception("[I] The computation was interrupted by user");
}
}
//' 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)
......@@ -132,28 +212,28 @@ void computeQuantilesCpp(const NumericMatrix swe,
if (check_interrupt()) {
interrupt = true;
}
std::vector<std::vector<double>> out_cpp(num_pixels);
const std::vector<size_t>& time_indices = timePools_cpp[doy-1];
const size_t num_pooled_timestep = time_indices.size();
for(size_t pixel=0; pixel<num_pixels; pixel++) {
const std::vector<size_t>& pixels_indices = spatialPools_cpp[pixel];
const size_t num_pooled_pixels = pixels_indices.size();
std::vector<double> quantiles_data(num_pooled_pixels*num_pooled_timestep);
for (int i=0; i < num_pooled_pixels; i++) {
for (int j=0; j < num_pooled_timestep; j++) {
for (size_t i=0; i < num_pooled_pixels; i++) {
for (size_t j=0; j < num_pooled_timestep; j++) {
quantiles_data[i+j*num_pooled_pixels] = swe_cpp[(pixels_indices[i]-1) +
(time_indices[j]-1)*num_pixels];
(time_indices[j]-1)*num_pixels];
}
}
remove_nan(quantiles_data);
const std::vector<double> quantiles = quantile_cpp(quantiles_data,probs_cpp);
out_cpp[pixel] = quantiles;
}
std::ofstream outFile(outputPath+"/quantiles_"+std::to_string(doy)+".txt");
for (const auto &line : out_cpp) {
for (const auto &e : line) {
......
......@@ -40,7 +40,7 @@ public:
void ma(std::vector<double>& vect, const int ma_value) {
const std::vector<double> vect_copy{vect};
const size_t size = vect.size();
const int size = vect.size();
const int lag = (ma_value/2);
for(int i = 0; i < size; ++i) {
for(int l = -lag; l <= lag; ++l) {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment