WSL/SLF GitLab Repository

quantile_application.R 5.21 KB
Newer Older
Adrien Michel's avatar
Adrien Michel committed
1
2
3
4
5
6
7
8
9
10
apply.qm.internal <- function(data, doy.selection, source.path, target.path,
                              type, ncores)
{
  if(ncores > 1){
    doParallel::registerDoParallel(ncores)
    `%do%` <- foreach::`%dopar%`
  } else{
    `%do%` <- foreach::`%do%`
  }
  data.corrected <- data
11
  pixels <- c(1:length(data$data[,1]))
Adrien Michel's avatar
Adrien Michel committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
  for(doy in names(doy.selection)) {
    doy.subset <- doy.selection[[doy]]
    PDF.source <- as.matrix(data.table::fread(paste0(source.path,
                                                    "/quantiles_",doy,".txt"),
                                              header=F))
    PDF.dest <- as.matrix(data.table::fread(paste0(target.path,
                                                  "/quantiles_",doy,".txt"),
                                            header=F))
    if(length(PDF.source[,1]) > length(pixels)) {
      stop(paste0("File ",source.path, "/quantiles_",
                  doy, ".txt has the wrong number of pixels"))
    }
    if(length(PDF.dest[,1]) > length(pixels)) {
      stop(paste0("File ",target.path, "/quantiles_",
                  doy, ".txt has the wrong number of pixels"))
    }
    correction = foreach::foreach(p = pixels) %do% {
      if(length(PDF.source[p,]) != length(PDF.dest[p,])){
        stop("Source and target quantile files have different number of qunatiles")
      }
      corrected <- rep(NA,length(doy.subset))
      i <- 1
      for(d in doy.subset) {
        quantile <- which(PDF.source[p,] > data$data[p,d])[1]-1
36
        quantile <- bound.quantile(quantile, length(PDF.source[p,]))
Adrien Michel's avatar
Adrien Michel committed
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
        corrected[i] <- data$data[p,d] - PDF.source[p,quantile] + PDF.dest[p,quantile]
        i <- i+1
      }
      as.vector(corrected)
    }
    for(p in pixels) {
      i=1
      for(s in sel){
        data.corrected$data[p,s]=correction[[p]][i]
        i=i+1
      }
    }
  }
  if(ncores > 1){
    doParallel::stopImplicitCluster()
  }
  return(data.corrected)
}


#' Apply the trained QM model to a given dataset
#'
#' \code{apply.model} applies a trained QM model (see
#' \code{\link{compute.quantiles()}}) to the dataset given as input.
#'
62
63
64
65
66
67
68
69
#' @param in.data.file Dataset on which the QM model should be applied. Dataset
#' should be a RDS file in a vectdata format (see \code{\link{vectorize.list()}})
#' @param in.data.file Path of the file to write the corrected dataset. Dataset
#' will be a RDS file in a vectdata format (see \code{\link{vectorize.list()}})
#' @param source.path Directory containing the directories with quantile distributions of the
#' training and target datasets (see \code{\link{vectorize.list()}})
#' @param type R or cpp to choose between the R or the C++ implementation
#' (see \code{\link{apply.qm.internal()}} and \code{\link{applyQMCpp()}}).
Adrien Michel's avatar
Adrien Michel committed
70
71
72
73
74
75
#' C++ implementation is faster and recommended, given the package has been
#' installed with OpenMP support.
#' @param ncores number of cores to be used
#'
#' @author Adrien Michel, 2021 (WSL/SLF)
#' @export
76
apply.model <- function(in.data.file, out.data.file, quantile.path, type, ncores)
Adrien Michel's avatar
Adrien Michel committed
77
{
78
79
  print("[I] Reading data")
  data <- readRDS(in.data.file)
Adrien Michel's avatar
Adrien Michel committed
80

81
  print("[I] Computing days of the year")
Adrien Michel's avatar
Adrien Michel committed
82
83
84
85
86
87
88
89
  doys=as.integer(strftime(data$time.date, format = "%j"))
  doys.unique=sort(unique(doys))
  `%do%` <- foreach::`%do%`
  doy.selection = foreach::foreach(doy=doys.unique) %do% {
    which(doys == doy)
  }
  names(doy.selection)=doys.unique

90
91
  print("[I] Correcting data")
  start.time <- Sys.time()
Adrien Michel's avatar
Adrien Michel committed
92
93
94
  if(type == "cpp") {
    corrected.swe <- applyQMCpp(data$data,
                                doySelection = doy.selection,
95
96
                                sourcePath = paste0(quantile.path,"/training_quantiles"),
                                targetPath = paste0(quantile.path,"/target_quantiles"),
Adrien Michel's avatar
Adrien Michel committed
97
98
99
100
101
102
                                ncores = ncores)
    data.corrected <- data
    data.corrected$data <- matrix(corrected.swe, ncol=dim(data$data)[2])
  } else if (type == "R"){
      data.corrected <- apply.qm.internal(data = data,
                                          doy.selection = doy.selection,
103
104
                                          source.path = paste0(quantile.path,"/training_quantiles"),
                                          target.path = paste0(quantile.path,"/target_quantiles"),
Adrien Michel's avatar
Adrien Michel committed
105
106
107
108
                                          ncores = ncores)
  } else {
    stop("Type sould be either cpp or R")
  }
109
110
111
112
113
114
115
116
117
  print(paste("[I] Correction computed in",
              as.numeric(round(difftime(Sys.time(),start.time, units = "mins"),2)),
              "minutes"))

  saveRDS(data.corrected,out.data.file)

  print("[I] Cleaning memory")
  rm(list=ls())
  invisible(gc())
Adrien Michel's avatar
Adrien Michel committed
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
}


#' Bound a given quantile
#'
#' \code{bound.quantile} bounds the input quantile between 2 and length of the
#' quantile distribution -1. Input value of NA are set to the length of the
#' quantile distribution -1 (see \code{\link{apply.qm.internal()}}, where values
#' bigger than any value in the quantile distribution are set to NA).
#'
#' @param quantile Quantile to be bounded
#' @param size Size of the quantile distribution used
#'
#' @return Bounded quantile
#'
#' @author Adrien Michel, 2021 (WSL/SLF)
bound.quantile <- function(quantile, size)
{
  if(is.na(quantile)) {quantile=size-1}
  else if(quantile<2) {quantile=2}
  else if(quantile>size-1) {quantile=size-1}
  return(quantile)
}