From 37bff7ce511dc8e5d2123e45c2393ade8fc60cfb Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Tue, 14 May 2019 14:04:59 +0200 Subject: [PATCH 01/19] New function to compute the properties of the underlying attractor. Dynamical conditions for the bias correction. --- R/ProxiesAttractor.R | 87 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) create mode 100644 R/ProxiesAttractor.R diff --git a/R/ProxiesAttractor.R b/R/ProxiesAttractor.R new file mode 100644 index 00000000..73f45227 --- /dev/null +++ b/R/ProxiesAttractor.R @@ -0,0 +1,87 @@ +#'ProxiesAttractor +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@description This function computes two dinamical proxies of the attractor: The local dimension (d) and the inverse of the persistence (theta). +#'These two parameters will be used as a condition for the bias correction function DynBiasCorrection +#'Funtion based on the matlab code (@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}) used in "Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#'The hammam effect or how a warm ocean enhances large scale atmospheric predictability. Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#' +#'@param dat data to create the attractor. Must be a matrix with the timesteps in rows and gridpoints in columns -> dat(time,grids) +# +#'@param quanti threshold for the computation of the attractor +#' +#'@param iplot to have a quick view of the proxies and the attractor +#' +#'@return dim and theta +#' +#'@return A list of 2: +#'\itemize{ +#'\item\code{$dim} {A vector with the same length than the timesteps: nrow(dat)} +#'\item\code{$theta} {A vector with the same length than the timesteps: nrow(dat)}} +#'@references +#' \url{https://doi.org/10.1038/s41467-019-09305-8} +#' +#'@import s2dverification +#'@import fields +#'@import pracma +#' +#'@examples +#' Creating an example of matrix dat(time,grids): +tm=24*6 # time +gm=5000 # grids +m=matrix(rand(1,tm*gm),nrow=tm,ncol=gm) +qm=0.90 # imposing a threshold +#'# selecting the option to create a plot to see the attractor 'iplot=TRUE' +#' example +ExampleAttractor=ProxiesAttractor(dat=m,quanti=qm,iplot=TRUE) +#' +ProxiesAttractor <- function(dat,quanti,iplot=FALSE){ + npoints=nrow(dat) + dim=numeric(npoints) + theta=numeric(npoints) + time=seq(1:npoints) + l=1 + for (l in 1:npoints){ + #computation of the distances from the iteration l + distance=pdist2(dat[l,],dat) + #Taking the logarithm of the distance + logdista=-log(distance) + n=1 + for(n in 1:npoints){ + if(logdista[n]=="Inf"){logdista[n]=NaN} + n=n+1} + #Compute the thheshold corresponding to the quantile + thresh=quantile(logdista, quanti,na.rm=TRUE) + + #Computation of the inverse of persistence 'theta' + Li<-which(logdista>thresh) + #Length of each cluster + Ti<-diff(Li) + N=length(Ti) + q=1-quanti + Si=Ti-1 + Nc=length(which(Si>0)) + N=length(Ti) + theta[l]=(sum(q*Si)+N+Nc-sqrt(((sum(q*Si)+N+Nc)^2)-8*Nc*sum(q*Si)))/(2*sum(q*Si)) + + #computation of local dimension 'dim' + #Sort the exceedances + logdista=sort(logdista) + #Find all the Peaks over Thresholds. + findidx=which(logdista>thresh) + logextr=logdista[findidx[[1]]:(length(logdista)-1)] + #The inverse of the dimension is just the average of the exceedances + dim[l]=1/mean(logextr-thresh) + l=l+1 + if(l Date: Mon, 17 May 2021 13:45:20 +0000 Subject: [PATCH 02/19] Update ProxiesAttractor.R --- R/ProxiesAttractor.R | 192 ++++++++++++++++++++++++++++--------------- 1 file changed, 126 insertions(+), 66 deletions(-) diff --git a/R/ProxiesAttractor.R b/R/ProxiesAttractor.R index 73f45227..f502cec1 100644 --- a/R/ProxiesAttractor.R +++ b/R/ProxiesAttractor.R @@ -1,87 +1,147 @@ #'ProxiesAttractor #' #'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} -#'@description This function computes two dinamical proxies of the attractor: The local dimension (d) and the inverse of the persistence (theta). -#'These two parameters will be used as a condition for the bias correction function DynBiasCorrection -#'Funtion based on the matlab code (@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr}) used in "Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). -#'The hammam effect or how a warm ocean enhances large scale atmospheric predictability. Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} #' -#'@param dat data to create the attractor. Must be a matrix with the timesteps in rows and gridpoints in columns -> dat(time,grids) -# -#'@param quanti threshold for the computation of the attractor +#'@description This function computes two dinamical proxies of the attractor: +#'The local dimension (d) and the inverse of the persistence (theta). +#'These two parameters will be used as a condition for the computation of dynamical +#'scores to measure predictability and to compute bias correction conditioned by +#'the dynamics with the function DynBiasCorrection +#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in +#'"Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#' Davide Faranda, Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. #' -#'@param iplot to have a quick view of the proxies and the attractor +#'@param data data to create the attractor. Must be a matrix with the timesteps in nrow +#'and the grids in ncol(dat(time,grids) +# +#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to +#' be provided as a list of length 2 with longitudes and latitudes +#' +#'@param iplot FALSE by default. If TRUE the function returns a simple plot +#'with three pannels on the top left a plot for local dimension 'd', on the top +#'right a plot for the inverse of the persistence 'theta', on the bottom a plot +#'for the 'attractor' plotted with the two properties 'd' (x axis) and 'theta' +#'(y axis) #' #'@return dim and theta -#' -#'@return A list of 2: -#'\itemize{ -#'\item\code{$dim} {A vector with the same length than the timesteps: nrow(dat)} -#'\item\code{$theta} {A vector with the same length than the timesteps: nrow(dat)}} -#'@references -#' \url{https://doi.org/10.1038/s41467-019-09305-8} -#' #'@import s2dverification #'@import fields -#'@import pracma #' #'@examples -#' Creating an example of matrix dat(time,grids): -tm=24*6 # time -gm=5000 # grids -m=matrix(rand(1,tm*gm),nrow=tm,ncol=gm) -qm=0.90 # imposing a threshold -#'# selecting the option to create a plot to see the attractor 'iplot=TRUE' -#' example -ExampleAttractor=ProxiesAttractor(dat=m,quanti=qm,iplot=TRUE) +#'# Example 1: computing the attractor using simple data +#'# Creating an example of matrix data(time,grids): +#'tm=2*6*3 # time +#'gm=4*10 # grids +#'m=matrix(rand(1,tm*gm),nrow=tm,ncol=gm) +#'qm=0.90 # imposing a threshold +#'Attractor=ProxiesAttractor(data=m,quanti=qm,iplot=FALSE) +#' +#'# Example 2: computing the attractor using simple s2dv data +#'expL <- rnorm(1:200) +#'dim(expL) <- c(member=10,lat = 4, lon = 5) +#'obsL <- c(rnorm(1:180),expL[1,,]*1.2) +#'dim(obsL) <- c(time = 10,lat = 4, lon = 5) +#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-") +#'time_expL <- time_obsL[1] +#'lon <- seq(-1,5,1.5) +#'lat <- seq(30,35,1.5) +#'qm=0.98 # too high for this short dataset, it is possible that doesn't +#'get the requirement, in that case it would be necessary select a lower qm +#'for instance qm=0.60 +#'qm=0.60 +#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' Dates = list(start = time_expL, end = time_expL)) +#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' Dates = list(start = time_obsL, end = time_obsL)) +#' +#'attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm) +#' #' -ProxiesAttractor <- function(dat,quanti,iplot=FALSE){ - npoints=nrow(dat) - dim=numeric(npoints) - theta=numeric(npoints) - time=seq(1:npoints) - l=1 - for (l in 1:npoints){ - #computation of the distances from the iteration l - distance=pdist2(dat[l,],dat) - #Taking the logarithm of the distance - logdista=-log(distance) - n=1 - for(n in 1:npoints){ - if(logdista[n]=="Inf"){logdista[n]=NaN} - n=n+1} +#' +#' +CST_ProxiesAttractor <- function(data,quanti,ncores=NULL){ + + if (!inherits(data, 's2dv_cube')) { + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (any(names(dim(data$data)) %in% 'sdate')) { + if (any(names(dim(data$data)) %in% 'ftime')) { + data$data <- MergeDims(data$data, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(data$data)) %in% 'time'))){ + stop('dimension time is mandatory') + } + if (any(names(dim(data$data)) %in% 'lat')) { + if (any(names(dim(data$data)) %in% 'lon')) { + data$data <- MergeDims(data$data, merge_dims = c('lon', 'lat'), + rename_dim = 'grid') + } + } + if(!(any(names(dim(data$data)) %in% 'grid'))){ + stop('dimension grid is mandatory') + } + + attractor <-Apply(data$data,target_dims=c('time','grid'), fun=ProxiesAttractor, + quanti=quanti,iplot=FALSE, ncores=ncores) + # rename dimensions + attractor <- lapply(attractor, FUN=function(x,dimname){names(dim(x))[dimname] <- 'time' + return(x)}, dimname = which(names(dim(attractor[[1]]))=='dim2')) + + return(list(dim=attractor$dim,theta=attractor$theta)) +} + +ProxiesAttractor <- function(data,quanti, iplot=FALSE,ncores=NULL){ + + logdista <- apply(data,1,function(x,y){ + -log(colMeans((y-x)^2))}, + y=t(data)) + + #Computation of theta + Theta <- function(logdista,quanti){ #Compute the thheshold corresponding to the quantile - thresh=quantile(logdista, quanti,na.rm=TRUE) - - #Computation of the inverse of persistence 'theta' - Li<-which(logdista>thresh) + thresh <- quantile(logdista,quanti,na.rm=TRUE) + logdista[which(logdista=='Inf')] <- NaN + Li <- which(as.vector(logdista) > as.numeric(thresh)) #Length of each cluster - Ti<-diff(Li) - N=length(Ti) - q=1-quanti - Si=Ti-1 - Nc=length(which(Si>0)) - N=length(Ti) - theta[l]=(sum(q*Si)+N+Nc-sqrt(((sum(q*Si)+N+Nc)^2)-8*Nc*sum(q*Si)))/(2*sum(q*Si)) - - #computation of local dimension 'dim' + Ti <- diff(Li) + N <- length(Ti) + q <- 1-quanti + Si <- Ti-1 + Nc <- length(which(Si>0)) + N <- length(Ti) + theta <- (sum(q*Si)+N+Nc-sqrt(((sum(q*Si)+N+Nc)^2)-8*Nc*sum(q*Si)))/(2*sum(q*Si)) #Sort the exceedances - logdista=sort(logdista) - #Find all the Peaks over Thresholds. - findidx=which(logdista>thresh) - logextr=logdista[findidx[[1]]:(length(logdista)-1)] + logdista <- sort(logdista) + #Find all the Peaks over Thresholds. + findidx <- which(as.vector(logdista) > as.numeric(thresh)) + if(length(findidx) < 1){ + stop("quanti too high for the length of the dataset") + } + logextr <- logdista[findidx[[1]]:(length(logdista)-1)] #The inverse of the dimension is just the average of the exceedances - dim[l]=1/mean(logextr-thresh) - l=l+1 - if(l Date: Tue, 18 May 2021 07:54:45 +0000 Subject: [PATCH 03/19] Update ProxiesAttractor.R --- R/ProxiesAttractor.R | 90 ++++++++++++-------------------------------- 1 file changed, 24 insertions(+), 66 deletions(-) diff --git a/R/ProxiesAttractor.R b/R/ProxiesAttractor.R index f502cec1..c400b98c 100644 --- a/R/ProxiesAttractor.R +++ b/R/ProxiesAttractor.R @@ -1,5 +1,7 @@ -#'ProxiesAttractor + +#'@rdname ProxiesAttractor #' +#'@title Computing two dinamical proxies of the attractor. #'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} #'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} #'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} @@ -9,12 +11,12 @@ #'The local dimension (d) and the inverse of the persistence (theta). #'These two parameters will be used as a condition for the computation of dynamical #'scores to measure predictability and to compute bias correction conditioned by -#'the dynamics with the function DynBiasCorrection -#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in -#'"Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#'the dynamics with the function DynBiasCorrection. +#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). #' The hammam effect or how a warm ocean enhances large scale atmospheric predictability. #' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " -#' Davide Faranda, Gabriele Messori and Pascal Yiou. (2017). +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). #' Dynamical proxies of North Atlantic predictability and extremes. #' Scientific Reports, 7-41278, 2017. #' @@ -24,6 +26,8 @@ #'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to #' be provided as a list of length 2 with longitudes and latitudes #' +#'@param ncores The number of cores to use in parallel computation +#' #'@param iplot FALSE by default. If TRUE the function returns a simple plot #'with three pannels on the top left a plot for local dimension 'd', on the top #'right a plot for the inverse of the persistence 'theta', on the bottom a plot @@ -32,10 +36,11 @@ #' #'@return dim and theta #'@import s2dverification +#'@import multiApply #'@import fields #' #'@examples -#'# Example 1: computing the attractor using simple data +#'# Example 1: Computing the attractor using simple data #'# Creating an example of matrix data(time,grids): #'tm=2*6*3 # time #'gm=4*10 # grids @@ -43,68 +48,21 @@ #'qm=0.90 # imposing a threshold #'Attractor=ProxiesAttractor(data=m,quanti=qm,iplot=FALSE) #' -#'# Example 2: computing the attractor using simple s2dv data -#'expL <- rnorm(1:200) -#'dim(expL) <- c(member=10,lat = 4, lon = 5) -#'obsL <- c(rnorm(1:180),expL[1,,]*1.2) -#'dim(obsL) <- c(time = 10,lat = 4, lon = 5) -#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-") -#'time_expL <- time_obsL[1] -#'lon <- seq(-1,5,1.5) -#'lat <- seq(30,35,1.5) -#'qm=0.98 # too high for this short dataset, it is possible that doesn't -#'get the requirement, in that case it would be necessary select a lower qm -#'for instance qm=0.60 -#'qm=0.60 -#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, -#' Dates = list(start = time_expL, end = time_expL)) -#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, -#' Dates = list(start = time_obsL, end = time_obsL)) -#' -#'attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm) -#' -#' -#' -#' -CST_ProxiesAttractor <- function(data,quanti,ncores=NULL){ - - if (!inherits(data, 's2dv_cube')) { - stop("Parameter 'data' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.") - } - if (any(names(dim(data$data)) %in% 'sdate')) { - if (any(names(dim(data$data)) %in% 'ftime')) { - data$data <- MergeDims(data$data, merge_dims = c('ftime', 'sdate'), - rename_dim = 'time') - } - } - if (!(any(names(dim(data$data)) %in% 'time'))){ - stop('dimension time is mandatory') - } - if (any(names(dim(data$data)) %in% 'lat')) { - if (any(names(dim(data$data)) %in% 'lon')) { - data$data <- MergeDims(data$data, merge_dims = c('lon', 'lat'), - rename_dim = 'grid') - } - } - if(!(any(names(dim(data$data)) %in% 'grid'))){ - stop('dimension grid is mandatory') - } - - attractor <-Apply(data$data,target_dims=c('time','grid'), fun=ProxiesAttractor, - quanti=quanti,iplot=FALSE, ncores=ncores) - # rename dimensions - attractor <- lapply(attractor, FUN=function(x,dimname){names(dim(x))[dimname] <- 'time' - return(x)}, dimname = which(names(dim(attractor[[1]]))=='dim2')) - - return(list(dim=attractor$dim,theta=attractor$theta)) -} + ProxiesAttractor <- function(data,quanti, iplot=FALSE,ncores=NULL){ - - logdista <- apply(data,1,function(x,y){ - -log(colMeans((y-x)^2))}, - y=t(data)) + if (is.null(data)) { + stop("Parameter 'data' is mandatory") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' is mandatory") + } + if (is.null(iplot)) { + stop("Parameter 'iplot' is mandatory and it is set to FALSE by default") + } + logdista <- Apply(data,target_dims =2, + fun = function(x,y){-log(colMeans((y-as.vector(x))^2))}, + y=t(data),ncores=ncores)[[1]] #Computation of theta Theta <- function(logdista,quanti){ -- GitLab From 6d2ecf5a406d1d89e1f01cffad33d7ec1d3cf078 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Tue, 18 May 2021 07:55:35 +0000 Subject: [PATCH 04/19] Add new file CST_ProxiesAttractor --- R/CST_ProxiesAttractor | 94 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 R/CST_ProxiesAttractor diff --git a/R/CST_ProxiesAttractor b/R/CST_ProxiesAttractor new file mode 100644 index 00000000..5ec3373e --- /dev/null +++ b/R/CST_ProxiesAttractor @@ -0,0 +1,94 @@ +#'@rdname CST_ProxiesAttractor +#'@title Computing two dinamical proxies of the attractor in s2dvcube. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function computes two dinamical proxies of the attractor: +#'The local dimension (d) and the inverse of the persistence (theta) for an +#''s2dvcube' object. +#'These two parameters will be used as a condition for the computation of dynamical +#'scores to measure predictability and to compute bias correction conditioned by +#'the dynamics with the function DynBiasCorrection +#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param data data to create the attractor. Must be a matrix with the timesteps in nrow +#'and the grids in ncol(dat(time,grids) +# +#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to +#' be provided as a list of length 2 with longitudes and latitudes +#' +#'@param ncores The number of cores to use in parallel computation +#' +#'@return dim and theta +#'@import s2dverification +#'@import fields +#'@import multiApply +#' +#'@examples +#'# Example 1: Computing the attractor using simple s2dv data +#'expL <- rnorm(1:200) +#'dim(expL) <- c(member=10,lat = 4, lon = 5) +#'obsL <- c(rnorm(1:180),expL[1,,]*1.2) +#'dim(obsL) <- c(time = 10,lat = 4, lon = 5) +#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-") +#'time_expL <- time_obsL[1] +#'lon <- seq(-1,5,1.5) +#'lat <- seq(30,35,1.5) +#'qm=0.98 # too high for this short dataset, it is possible that doesn't +#'get the requirement, in that case it would be necessary select a lower qm +#'for instance qm=0.60 +#'qm=0.60 +#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' Dates = list(start = time_expL, end = time_expL)) +#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' Dates = list(start = time_obsL, end = time_obsL)) +#' +#'attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm) +#' +CST_ProxiesAttractor <- function(data,quanti,ncores=NULL){ + if (is.null(data)) { + stop("Parameter 'data' is mandatory") + } + if (!inherits(data, 's2dv_cube')) { + stop("Parameter 'data' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' is mandatory") + } + if (any(names(dim(data$data)) %in% 'sdate')) { + if (any(names(dim(data$data)) %in% 'ftime')) { + data$data <- MergeDims(data$data, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(data$data)) %in% 'time'))){ + stop('dimension time is mandatory') + } + if (any(names(dim(data$data)) %in% 'lat')) { + if (any(names(dim(data$data)) %in% 'lon')) { + data$data <- MergeDims(data$data, merge_dims = c('lon', 'lat'), + rename_dim = 'grid') + } + } + if(!(any(names(dim(data$data)) %in% 'grid'))){ + stop('dimension grid is mandatory') + } + + attractor <-Apply(data$data,target_dims=c('time','grid'), fun=ProxiesAttractor, + quanti=quanti,iplot=FALSE, ncores=ncores) + # rename dimensions + attractor <- lapply(attractor, FUN=function(x,dimname){names(dim(x))[dimname] <- 'time' + return(x)}, dimname = which(names(dim(attractor[[1]]))=='dim2')) + + return(list(dim=attractor$dim,theta=attractor$theta)) +} -- GitLab From c62db3be28119f2c99b0de7976d0876904ac9da2 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Tue, 18 May 2021 07:55:56 +0000 Subject: [PATCH 05/19] Update CST_ProxiesAttractor --- R/{CST_ProxiesAttractor => CST_ProxiesAttractor.R} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename R/{CST_ProxiesAttractor => CST_ProxiesAttractor.R} (100%) diff --git a/R/CST_ProxiesAttractor b/R/CST_ProxiesAttractor.R similarity index 100% rename from R/CST_ProxiesAttractor rename to R/CST_ProxiesAttractor.R -- GitLab From 41a5bb25e17246531717f6d356185df3b30a00cd Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Fri, 21 May 2021 16:56:03 +0000 Subject: [PATCH 06/19] Add new file --- R/CST_Predictability.R | 1 + 1 file changed, 1 insertion(+) create mode 100644 R/CST_Predictability.R diff --git a/R/CST_Predictability.R b/R/CST_Predictability.R new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/R/CST_Predictability.R @@ -0,0 +1 @@ + -- GitLab From 01a42b6717d6287527b05d4150a6ca278d7d9cd6 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Fri, 21 May 2021 16:56:33 +0000 Subject: [PATCH 07/19] Update CST_Predictability.R --- R/CST_Predictability.R | 118 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 118 insertions(+) diff --git a/R/CST_Predictability.R b/R/CST_Predictability.R index 8b137891..fffa1c64 100644 --- a/R/CST_Predictability.R +++ b/R/CST_Predictability.R @@ -1 +1,119 @@ +#'@rdname CST_Predictability +#'@title Computing scores of predictability using two dinamical proxies +#'based on dynamical systems theory. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function divides in terciles the two dynamical proxies +#'computed with CST_ProxiesAttractor or ProxiesAttractor. These terciles will +#'be used to measure the predictability of the system in dyn_scores. When the +#'local dimension 'dim' is small and the inverse of persistence 'theta' is +#'small the predictability is higher, and viceversa. +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param dim data to create the attractor. Must be a matrix with the timesteps in nrow and the grids in ncol +#'for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints +# +#'@param theta list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes +#' +#'@param ncores The number of cores to use in parallel computation +#' +#'@return pred.dim a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +#'contains values of local dimension 'dim' divided by terciles: +#'d1: lower tercile and more predictability, +#'d2: middle tercile, +#'d3: higher tercile and less predictability +#'The 'pos.d' list contains the position of each tercile in parameter 'dim' +#' +#'@return pred.theta a list of two lists 'qtheta' and 'pos.t'. +#'The 'qtheta' list contains values of the inverse of persistence 'theta' +#'divided by terciles: +#'th1: lower tercile and more predictability, +#'th2: middle tercile, +#'th3: higher tercile and less predictability +#'The 'pos.t' list contains the position of each tercile in parameter 'theta' +#' +#'@return dyn_scores values from 0 to 1. A dyn_score of 1 indicates higher +#'predictability. +#' +#'@import s2dverification +#'@import multiApply +#'@import CSTools +#' +#'@examples +#'# Creating an example of matrix dat(time,grids): +#' m=matrix(rand(1,2000)*10,nrow=50,ncol=40) +#'# imposing a threshold +#' quanti=0.90 +#'# computing dyn_scores from parameters dim and theta of the attractor +#' attractor = ProxiesAttractor(dat=m,quanti=0.60,iplot=FALSE) +#' predyn = CST_Predictability(dim=attractor$dim,theta=attractor$theta) +#' +CST_Predictability <- function(dim,theta,ncores=NULL){ + if (is.null(dim)) { + stop("Parameter 'dim' is mandatory") + } + if (is.null(theta)) { + stop("Parameter 'theta' is mandatory") + } + if (length(dim)!=length(theta)) { + stop("Parameters 'dim' and 'theta' must have the same length") + } + pos <- c(1:length(dim)) + + # dim + qd1 <- quantile(dim,0.33,na.rm=TRUE) + qd3 <- quantile(dim,0.66,na.rm=TRUE) + d3 <- which(dim>=qd3) + d1 <- which(dim<=qd1) + d2 <- which(dim>qd1 & dim=qt3) + th1 <- which(theta<=qt1) + th2 <- which(theta>qt1 & theta Date: Fri, 21 May 2021 17:25:13 +0000 Subject: [PATCH 08/19] Add new file --- R/DynBiasCorrection.R | 92 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) create mode 100644 R/DynBiasCorrection.R diff --git a/R/DynBiasCorrection.R b/R/DynBiasCorrection.R new file mode 100644 index 00000000..36f7864f --- /dev/null +++ b/R/DynBiasCorrection.R @@ -0,0 +1,92 @@ +#'@rdname DynBiasCorrection +#'@title Performing a Bias Correction conditioned by the dynamical +#'properties of the data. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#' +#'@description This function perform a bias correction conditioned by the +#'dynamical properties of the dataset. This function used the functions +#''CST_Predictability' to divide in terciles the two dynamical proxies +#'computed with CST_ProxiesAttractor or ProxiesAttractor. A bias correction +#'between the model and the observations is performed using the division into +#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +#'For instance, model values with lower 'dim' will be corrected with observed +#'values with lower 'dim', and the same for theta. The function gives two options +#'of bias correction: one for 'dim' and one for 'theta' +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param pred.dim.mod output of CST_Predictability or a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +#'contains values of local dimension 'dim' divided by terciles: +#'d1: lower tercile and more predictability, +#'d2: middle tercile, +#'d3: higher tercile and less predictability +#'The 'pos.d' list contains the position of each tercile in parameter 'dim' +#'@param pred.theta.mod output of CST_Predictability or a list of two lists 'qtheta' and 'pos.t'. +#'The 'qtheta' list contains values of the inverse of persistence 'theta' +#'divided by terciles: +#'th1: lower tercile and more predictability, +#'th2: middle tercile, +#'th3: higher tercile and less predictability +#'The 'pos.t' list contains the position of each tercile in parameter 'theta' +#' +#'@param pred.dim.obs same than pred.dim.obs but for the observations +#' +#'@param pred.theta.obs same than pred.theta.obs but for the observations +#' +#'@param ncores The number of cores to use in parallel computation +#' +#'@return dim.bias an s2dvcube object with a bias correction performed +#'conditioned by local dimension 'dim' +#' +#'@return theta.bias an s2dvcube object with a bias correction performed +#'conditioned by the inverse of the persistence 'theta' +#' +#'@import s2dverification +#'@import multiApply +#'@import CSTools +#' +#'@examples +#'# example 1: simple data s2dvcube style +# expL <- rnorm(1:200) +# dim(expL) <- c(member=10,lat = 4, lon = 5) +# obsL <- c(rnorm(1:180),expL[1,,]*1.2) +# dim(obsL) <- c(time = 10,lat = 4, lon = 5) +# time_obsL <- paste(rep("01", 50), rep("01", 50), 1953 : 2003, sep = "-") +# time_expL <- time_obsL[1] +# lon <- seq(-1,5,1.5) +# lat <- seq(30,35,1.5) +# qm=0.98 +# expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +# Dates = list(start = time_expL, end = time_expL)) +# obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +# Dates = list(start = time_obsL, end = time_obsL)) +# attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm) +# predynObs = CST_Predictability(dim=attractor$dim,theta=attractor$theta) +# attractorMod<- CST_ProxiesAttractor(data=expL,quanti=qm) +# predynMod = CST_Predictability(dim=attractorMod$dim,theta=attractorMod$theta) +# dynbias= DynBiasCorrection(pred.dim.obs=predynObs$pred.dim,pred.theta.obs=predynMod$pred.theta, +# pred.dim.mod=predynObs$pred.dim,pred.theta.mod=predynMod$pred.theta) + +DynBiasCorrection<- function(pred.dim.obs,pred.theta.obs,pred.dim.mod, + pred.theta.mod,ncores=NULL){ + if (is.null(pred.dim.mod)) { + stop("Parameter 'pred.dim.mod' is mandatory") + } + if (is.null(pred.dim.obs)) { + stop("Parameter 'pred.dim.obs' is mandatory") + } + if (is.null(pred.theta.mod)) { + stop("Parameter 'pred.dim.mod' is mandatory") + } + if (is.null(pred.theta.obs)) { + stop("Parameter 'pred.dim.obs' is mandatory") + } + return(list(dim.bias=dim.bias,theta.bias=theta.bias)) +} -- GitLab From b5693b0be59fac55bff02645ef988fb8854b37e1 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 31 May 2021 19:14:47 +0200 Subject: [PATCH 09/19] Changes to follow guidelines --- DESCRIPTION | 6 +- NAMESPACE | 6 + R/CST_ProxiesAttractor.R | 191 ++++++++++++++----- R/DynBiasCorrection.R | 41 ++-- R/{CST_Predictability.R => Predictability.R} | 70 ++++--- man/AnoMultiMetric.Rd | 3 +- man/CST_ProxiesAttractor.Rd | 52 +++++ man/DynBiasCorrection.Rd | 89 +++++++++ man/Predictability.Rd | 72 +++++++ man/ProxiesAttractor.Rd | 95 +++++++++ 10 files changed, 516 insertions(+), 109 deletions(-) rename R/{CST_Predictability.R => Predictability.R} (72%) create mode 100644 man/CST_ProxiesAttractor.Rd create mode 100644 man/DynBiasCorrection.Rd create mode 100644 man/Predictability.Rd create mode 100644 man/ProxiesAttractor.Rd diff --git a/DESCRIPTION b/DESCRIPTION index c7355509..3c22d1e5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -21,10 +21,12 @@ Depends: Imports: s2dverification, abind, - stats + stats, + multiApply, + fields License: LGPL-3 URL: https://earth.bsc.es/gitlab/external/cstools BugReports: https://earth.bsc.es/gitlab/external/cstools/issues Encoding: UTF-8 LazyData: true -RoxygenNote: 6.0.1.9000 +RoxygenNote: 7.0.2 diff --git a/NAMESPACE b/NAMESPACE index 576f3a5d..5bb49e11 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,12 @@ # Generated by roxygen2: do not edit by hand export(AnoMultiMetric) +export(CST_ProxiesAttractor) +export(DynBiasCorrection) export(MultivarRMSE) +export(Predictability) +export(ProxiesAttractor) +import(fields) +import(multiApply) import(s2dverification) import(stats) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 5ec3373e..903f05ca 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -1,5 +1,5 @@ #'@rdname CST_ProxiesAttractor -#'@title Computing two dinamical proxies of the attractor in s2dvcube. +#'@title Computing two dinamical proxies of the attractor in s2dv_cube. #' #'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} #'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} @@ -8,7 +8,7 @@ #' #'@description This function computes two dinamical proxies of the attractor: #'The local dimension (d) and the inverse of the persistence (theta) for an -#''s2dvcube' object. +#''s2dv_cube' object. #'These two parameters will be used as a condition for the computation of dynamical #'scores to measure predictability and to compute bias correction conditioned by #'the dynamics with the function DynBiasCorrection @@ -20,7 +20,7 @@ #' Dynamical proxies of North Atlantic predictability and extremes. #' Scientific Reports, 7-41278, 2017. #' -#'@param data data to create the attractor. Must be a matrix with the timesteps in nrow +#'@param data a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow #'and the grids in ncol(dat(time,grids) # #'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to @@ -29,66 +29,165 @@ #'@param ncores The number of cores to use in parallel computation #' #'@return dim and theta -#'@import s2dverification -#'@import fields -#'@import multiApply #' #'@examples #'# Example 1: Computing the attractor using simple s2dv data -#'expL <- rnorm(1:200) -#'dim(expL) <- c(member=10,lat = 4, lon = 5) -#'obsL <- c(rnorm(1:180),expL[1,,]*1.2) -#'dim(obsL) <- c(time = 10,lat = 4, lon = 5) -#'time_obsL <- paste(rep("01", 10), rep("01", 10), 1994 : 2003, sep = "-") -#'time_expL <- time_obsL[1] -#'lon <- seq(-1,5,1.5) -#'lat <- seq(30,35,1.5) -#'qm=0.98 # too high for this short dataset, it is possible that doesn't -#'get the requirement, in that case it would be necessary select a lower qm -#'for instance qm=0.60 -#'qm=0.60 -#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, -#' Dates = list(start = time_expL, end = time_expL)) -#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, -#' Dates = list(start = time_obsL, end = time_obsL)) -#' -#'attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm) +#'attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6) #' -CST_ProxiesAttractor <- function(data,quanti,ncores=NULL){ - if (is.null(data)) { - stop("Parameter 'data' is mandatory") - } +#'@export +CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ if (!inherits(data, 's2dv_cube')) { stop("Parameter 'data' must be of the class 's2dv_cube', ", "as output by CSTools::CST_Load.") } if (is.null(quanti)) { - stop("Parameter 'quanti' is mandatory") + stop("Parameter 'quanti' cannot be NULL.") } - if (any(names(dim(data$data)) %in% 'sdate')) { - if (any(names(dim(data$data)) %in% 'ftime')) { - data$data <- MergeDims(data$data, merge_dims = c('ftime', 'sdate'), + + data$data <- ProxiesAttractor(data = data$data, quanti = quanti, ncores = ncores) + + return(data) +} +#'@rdname ProxiesAttractor +#' +#'@title Computing two dinamical proxies of the attractor. +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function computes two dinamical proxies of the attractor: +#'The local dimension (d) and the inverse of the persistence (theta). +#'These two parameters will be used as a condition for the computation of dynamical +#'scores to measure predictability and to compute bias correction conditioned by +#'the dynamics with the function DynBiasCorrection. +#'Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +#' The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +#' Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param data a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'. +#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to +#' be provided as a list of length 2 with longitudes and latitudes. +#'@param iplot a logical indicating if to plot results. +#'@param ncores The number of cores to use in parallel computation. +#' +#'@param iplot FALSE by default. If TRUE the function returns a simple plot +#'with three pannels on the top left a plot for local dimension 'd', on the top +#'right a plot for the inverse of the persistence 'theta', on the bottom a plot +#'for the 'attractor' plotted with the two properties 'd' (x axis) and 'theta' +#'(y axis) +#' +#'@return dim and theta +#'@import multiApply +#' +#'@examples +#'# Example 1: Computing the attractor using simple data +#'# Creating an example of matrix data(time,grids): +#'mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) +#'qm <- 0.90 # imposing a threshold +#'Attractor <- ProxiesAttractor(data = mat, quanti = qm) +#' +#'@export + +ProxiesAttractor <- function(data, quanti, iplot = FALSE, ncores = NULL){ + if (is.null(data)) { + stop("Parameter 'data' cannot be NULL.") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' is mandatory") + } + if (!is.logical(iplot) || is.null(iplot)) { + stop("Parameter 'iplot' is required and needs to be TRUE or FALSE.") + } + if (any(names(dim(data)) %in% 'sdate')) { + if (any(names(dim(data)) %in% 'ftime')) { + data <- MergeDims(data, merge_dims = c('ftime', 'sdate'), rename_dim = 'time') } } - if (!(any(names(dim(data$data)) %in% 'time'))){ - stop('dimension time is mandatory') + if (!(any(names(dim(data)) %in% 'time'))){ + stop("Parameter 'data' must have a temporal dimension named 'time'.") } - if (any(names(dim(data$data)) %in% 'lat')) { - if (any(names(dim(data$data)) %in% 'lon')) { - data$data <- MergeDims(data$data, merge_dims = c('lon', 'lat'), + if (any(names(dim(data)) %in% 'lat')) { + if (any(names(dim(data)) %in% 'lon')) { + data <- MergeDims(data, merge_dims = c('lon', 'lat'), rename_dim = 'grid') } } - if(!(any(names(dim(data$data)) %in% 'grid'))){ - stop('dimension grid is mandatory') - } - - attractor <-Apply(data$data,target_dims=c('time','grid'), fun=ProxiesAttractor, - quanti=quanti,iplot=FALSE, ncores=ncores) + if (any(names(dim(data)) %in% 'latitude')) { + if (any(names(dim(data)) %in% 'longitude')) { + data <- MergeDims(data, merge_dims = c('longitude', 'latitude'), + rename_dim = 'grid') + } + } + if(!(any(names(dim(data)) %in% 'grid'))){ + stop("Parameter 'data' must have a spatial dimension named 'grid'.") + } + attractor <- Apply(data, target_dims = c('time', 'grid'), + fun = .proxiesattractor, + quanti = quanti , iplot = FALSE, ncores = ncores) # rename dimensions - attractor <- lapply(attractor, FUN=function(x,dimname){names(dim(x))[dimname] <- 'time' - return(x)}, dimname = which(names(dim(attractor[[1]]))=='dim2')) + attractor <- lapply(attractor, + FUN = function(x, dimname){ + names(dim(x))[dimname] <- 'time' + return(x)}, + dimname = which(names(dim(attractor[[1]])) == 'dim2')) + return(list(dim = attractor$dim, theta = attractor$theta)) +} + +.proxiesattractor <- function(data, quanti, iplot = FALSE) { + # expected dimensions data: time and grid + logdista <- Apply(data, target_dims = 'grid', + fun = function(x, y){ + -log(colMeans((y - as.vector(x))^2))}, + y = t(data))[[1]] + + #Computation of theta + Theta <- function(logdista, quanti){ + #Compute the thheshold corresponding to the quantile + thresh <- quantile(logdista, quanti, na.rm = TRUE) + logdista[which(logdista == 'Inf')] <- NaN + Li <- which(as.vector(logdista) > as.numeric(thresh)) + #Length of each cluster + Ti <- diff(Li) + N <- length(Ti) + q <- 1 - quanti + Si <- Ti - 1 + Nc <- length(which(Si > 0)) + N <- length(Ti) + theta <- (sum(q * Si) + N + Nc - sqrt(((sum(q * Si) + N + Nc)^2) - + 8 * Nc * sum(q * Si))) / (2 * sum(q * Si)) + #Sort the exceedances + logdista <- sort(logdista) + #Find all the Peaks over Thresholds. + findidx <- which(as.vector(logdista) > as.numeric(thresh)) + if(length(findidx) < 1) { + stop("Parameter 'quanti' is too high for the length of the data provided.") + } + logextr <- logdista[findidx[[1]]:(length(logdista) - 1)] + #The inverse of the dimension is just the average of the exceedances + dim <- 1 /mean(as.numeric(logextr) - as.numeric(thresh)) + return(list(dim = dim, theta = theta)) + } + names(dim(logdista)) <- c('dim1', 'dim2') + proxies <- Apply(data = list(logdista = logdista), + target_dims = list('dim1'), fun = Theta, quanti = quanti, + ncores = ncores) + if(iplot == TRUE) { + time = c(1:length(proxies$theta)) + layout(matrix(c(1, 3, 2, 3), 2, 2)) + plot(time, proxies$dim, xlab = 'time', ylab = 'd', + main = 'local dimension', type = 'l') + plot(time, proxies$theta, xlab = 'time', ylab = 'theta', main = 'theta') + plot(proxies$dim, proxies$theta, col = 'blue', + main = "Proxies of the Attractor", + xlab = "local dimension", ylab = "theta", lwd = 8, 'p') + } - return(list(dim=attractor$dim,theta=attractor$theta)) + return(list(dim = proxies$dim, theta = proxies$theta)) } + diff --git a/R/DynBiasCorrection.R b/R/DynBiasCorrection.R index 36f7864f..ccab645f 100644 --- a/R/DynBiasCorrection.R +++ b/R/DynBiasCorrection.R @@ -48,32 +48,29 @@ #'@return theta.bias an s2dvcube object with a bias correction performed #'conditioned by the inverse of the persistence 'theta' #' -#'@import s2dverification -#'@import multiApply -#'@import CSTools #' #'@examples #'# example 1: simple data s2dvcube style -# expL <- rnorm(1:200) -# dim(expL) <- c(member=10,lat = 4, lon = 5) -# obsL <- c(rnorm(1:180),expL[1,,]*1.2) -# dim(obsL) <- c(time = 10,lat = 4, lon = 5) -# time_obsL <- paste(rep("01", 50), rep("01", 50), 1953 : 2003, sep = "-") -# time_expL <- time_obsL[1] -# lon <- seq(-1,5,1.5) -# lat <- seq(30,35,1.5) -# qm=0.98 -# expL <- s2dv_cube(data = expL, lat = lat, lon = lon, -# Dates = list(start = time_expL, end = time_expL)) -# obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, -# Dates = list(start = time_obsL, end = time_obsL)) -# attractor<- CST_ProxiesAttractor(data=obsL,quanti=qm) -# predynObs = CST_Predictability(dim=attractor$dim,theta=attractor$theta) -# attractorMod<- CST_ProxiesAttractor(data=expL,quanti=qm) -# predynMod = CST_Predictability(dim=attractorMod$dim,theta=attractorMod$theta) -# dynbias= DynBiasCorrection(pred.dim.obs=predynObs$pred.dim,pred.theta.obs=predynMod$pred.theta, +#' expL <- rnorm(1:200) +#' dim(expL) <- c(member=10,lat = 4, lon = 5) +#' obsL <- c(rnorm(1:180),expL[1,,]*1.2) +#' dim(obsL) <- c(time = 10,lat = 4, lon = 5) +#' time_obsL <- paste(rep("01", 50), rep("01", 50), 1953 : 2003, sep = "-") +#' time_expL <- time_obsL[1] +#' lon <- seq(-1,5,1.5) +#' lat <- seq(30,35,1.5) +#' qm=0.98 +#' expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' Dates = list(start = time_expL, end = time_expL)) +#' obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' Dates = list(start = time_obsL, end = time_obsL)) +#' attractor <- CST_ProxiesAttractor(data=obsL,quanti=qm) +#' predynObs = Predictability(dim=attractor$dim,theta=attractor$theta) +#' attractorMod<- CST_ProxiesAttractor(data=expL,quanti=qm) +#' predynMod = CST_Predictability(dim=attractorMod$dim,theta=attractorMod$theta) +#' dynbias= DynBiasCorrection(pred.dim.obs=predynObs$pred.dim,pred.theta.obs=predynMod$pred.theta, # pred.dim.mod=predynObs$pred.dim,pred.theta.mod=predynMod$pred.theta) - +#'@export DynBiasCorrection<- function(pred.dim.obs,pred.theta.obs,pred.dim.mod, pred.theta.mod,ncores=NULL){ if (is.null(pred.dim.mod)) { diff --git a/R/CST_Predictability.R b/R/Predictability.R similarity index 72% rename from R/CST_Predictability.R rename to R/Predictability.R index fffa1c64..1989ef3d 100644 --- a/R/CST_Predictability.R +++ b/R/Predictability.R @@ -1,4 +1,4 @@ -#'@rdname CST_Predictability +#'@rdname Predictability #'@title Computing scores of predictability using two dinamical proxies #'based on dynamical systems theory. #' @@ -26,68 +26,64 @@ # #'@param theta list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes #' -#'@param ncores The number of cores to use in parallel computation -#' -#'@return pred.dim a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +#'@return A list of length 2: +#' \itemize{ +#' \item\code{pred.dim} {a list of two lists 'qdim' and 'pos.d'. The 'qdim' list #'contains values of local dimension 'dim' divided by terciles: #'d1: lower tercile and more predictability, #'d2: middle tercile, #'d3: higher tercile and less predictability -#'The 'pos.d' list contains the position of each tercile in parameter 'dim' +#'The 'pos.d' list contains the position of each tercile in parameter 'dim'} #' -#'@return pred.theta a list of two lists 'qtheta' and 'pos.t'. +#' \item\code{pred.theta} {a list of two lists 'qtheta' and 'pos.t'. #'The 'qtheta' list contains values of the inverse of persistence 'theta' #'divided by terciles: #'th1: lower tercile and more predictability, #'th2: middle tercile, #'th3: higher tercile and less predictability -#'The 'pos.t' list contains the position of each tercile in parameter 'theta' -#' +#'The 'pos.t' list contains the position of each tercile in parameter 'theta'} +#'} #'@return dyn_scores values from 0 to 1. A dyn_score of 1 indicates higher #'predictability. #' -#'@import s2dverification -#'@import multiApply -#'@import CSTools -#' #'@examples #'# Creating an example of matrix dat(time,grids): -#' m=matrix(rand(1,2000)*10,nrow=50,ncol=40) +#'m <- matrix(rnorm(2000) * 10, nrow = 50, ncol = 40) #'# imposing a threshold -#' quanti=0.90 +#' quanti <- 0.90 #'# computing dyn_scores from parameters dim and theta of the attractor -#' attractor = ProxiesAttractor(dat=m,quanti=0.60,iplot=FALSE) -#' predyn = CST_Predictability(dim=attractor$dim,theta=attractor$theta) -#' -CST_Predictability <- function(dim,theta,ncores=NULL){ +#' attractor <- ProxiesAttractor(dat = m, quanti = 0.60, iplot = FALSE) +#' predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) +#'@export +Predictability <- function(dim, theta) { if (is.null(dim)) { - stop("Parameter 'dim' is mandatory") + stop("Parameter 'dim' cannot be NULL.") } if (is.null(theta)) { - stop("Parameter 'theta' is mandatory") + stop("Parameter 'theta' cannot be NULL.") } - if (length(dim)!=length(theta)) { + if (length(dim) != length(theta)) { stop("Parameters 'dim' and 'theta' must have the same length") } pos <- c(1:length(dim)) # dim - qd1 <- quantile(dim,0.33,na.rm=TRUE) - qd3 <- quantile(dim,0.66,na.rm=TRUE) - d3 <- which(dim>=qd3) - d1 <- which(dim<=qd1) - d2 <- which(dim>qd1 & dim= qd3) + d1 <- which(dim <= qd1) + d2 <- which(dim > qd1 & dim < qd3) + qdim <- list(d1 = dim[d1], d2 = dim[d2], d3 = dim[d3]) pos.d <- list(d1=d1,d2=d2,d3=d3) #theta - qt1 <- quantile(theta,0.33,na.rm=TRUE) - qt3 <- quantile(theta,0.66,na.rm=TRUE) - th3 <- which(theta>=qt3) - th1 <- which(theta<=qt1) - th2 <- which(theta>qt1 & theta= qt3) + th1 <- which(theta <= qt1) + th2 <- which(theta > qt1 & theta < qt3) + qtheta <- list(th1 = theta[th1], th2 = theta[th2], th3 = theta[th3]) + pos.t <- list(th1 = th1, th2 = th2, th3 = th3) #scores position d1th1 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th1)] @@ -111,9 +107,9 @@ CST_Predictability <- function(dim,theta,ncores=NULL){ dyn_scores[which(pos %in% d3th2)]<- 2/9 dyn_scores[which(pos %in% d3th3)]<- 1/9 - return(list(pred.dim=list(qdim=qdim,pos.d=pos.d), - pred.theta=list(qtheta=qtheta,pos.t=pos.t), - dyn_scores=dyn_scores)) +return(list(pred.dimi =list(qdim = qdim, pos.d = pos.d), + pred.theta = list(qtheta = qtheta, pos.t = pos.t), + dyn_scores = dyn_scores)) } diff --git a/man/AnoMultiMetric.Rd b/man/AnoMultiMetric.Rd index 0071a721..28869b2c 100644 --- a/man/AnoMultiMetric.Rd +++ b/man/AnoMultiMetric.Rd @@ -4,8 +4,7 @@ \alias{AnoMultiMetric} \title{Multiple Metrics applied in Multiple Model Anomalies} \usage{ -AnoMultiMetric(data, metric = "correlation", multimodel = TRUE, - names = NULL) +AnoMultiMetric(data, metric = "correlation", multimodel = TRUE, names = NULL) } \arguments{ \item{data}{an s2dverification object (list) giving as output of \code{Load} function from S2dverification package.} diff --git a/man/CST_ProxiesAttractor.Rd b/man/CST_ProxiesAttractor.Rd new file mode 100644 index 00000000..9757ff10 --- /dev/null +++ b/man/CST_ProxiesAttractor.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_ProxiesAttractor.R +\name{CST_ProxiesAttractor} +\alias{CST_ProxiesAttractor} +\title{Computing two dinamical proxies of the attractor in s2dv_cube.} +\usage{ +CST_ProxiesAttractor(data, quanti, ncores = NULL) +} +\arguments{ +\item{data}{a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow +and the grids in ncol(dat(time,grids)} + +\item{quanti}{list of arbitrary length of secondary grids. Each secondary grid is to +be provided as a list of length 2 with longitudes and latitudes} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +dim and theta +} +\description{ +This function computes two dinamical proxies of the attractor: +The local dimension (d) and the inverse of the persistence (theta) for an +'s2dv_cube' object. +These two parameters will be used as a condition for the computation of dynamical +scores to measure predictability and to compute bias correction conditioned by +the dynamics with the function DynBiasCorrection +Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in +} +\examples{ +# Example 1: Computing the attractor using simple s2dv data +attractor <- CST_ProxiesAttractor(data = lonlat_data$obs, quanti = 0.6) + +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/DynBiasCorrection.Rd b/man/DynBiasCorrection.Rd new file mode 100644 index 00000000..7acfd1f0 --- /dev/null +++ b/man/DynBiasCorrection.Rd @@ -0,0 +1,89 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/DynBiasCorrection.R +\name{DynBiasCorrection} +\alias{DynBiasCorrection} +\title{Performing a Bias Correction conditioned by the dynamical +properties of the data.} +\usage{ +DynBiasCorrection( + pred.dim.obs, + pred.theta.obs, + pred.dim.mod, + pred.theta.mod, + ncores = NULL +) +} +\arguments{ +\item{pred.dim.obs}{same than pred.dim.obs but for the observations} + +\item{pred.theta.obs}{same than pred.theta.obs but for the observations} + +\item{pred.dim.mod}{output of CST_Predictability or a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +contains values of local dimension 'dim' divided by terciles: +d1: lower tercile and more predictability, +d2: middle tercile, +d3: higher tercile and less predictability +The 'pos.d' list contains the position of each tercile in parameter 'dim'} + +\item{pred.theta.mod}{output of CST_Predictability or a list of two lists 'qtheta' and 'pos.t'. +The 'qtheta' list contains values of the inverse of persistence 'theta' +divided by terciles: +th1: lower tercile and more predictability, +th2: middle tercile, +th3: higher tercile and less predictability +The 'pos.t' list contains the position of each tercile in parameter 'theta'} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +dim.bias an s2dvcube object with a bias correction performed +conditioned by local dimension 'dim' + +theta.bias an s2dvcube object with a bias correction performed +conditioned by the inverse of the persistence 'theta' +} +\description{ +This function perform a bias correction conditioned by the +dynamical properties of the dataset. This function used the functions +'CST_Predictability' to divide in terciles the two dynamical proxies +computed with CST_ProxiesAttractor or ProxiesAttractor. A bias correction +between the model and the observations is performed using the division into +terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +For instance, model values with lower 'dim' will be corrected with observed +values with lower 'dim', and the same for theta. The function gives two options +of bias correction: one for 'dim' and one for 'theta' +} +\examples{ +# example 1: simple data s2dvcube style +expL <- rnorm(1:200) +dim(expL) <- c(member=10,lat = 4, lon = 5) +obsL <- c(rnorm(1:180),expL[1,,]*1.2) +dim(obsL) <- c(time = 10,lat = 4, lon = 5) +time_obsL <- paste(rep("01", 50), rep("01", 50), 1953 : 2003, sep = "-") +time_expL <- time_obsL[1] +lon <- seq(-1,5,1.5) +lat <- seq(30,35,1.5) +qm=0.98 +expL <- s2dv_cube(data = expL, lat = lat, lon = lon, + Dates = list(start = time_expL, end = time_expL)) +obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, + Dates = list(start = time_obsL, end = time_obsL)) +attractor <- CST_ProxiesAttractor(data=obsL,quanti=qm) +predynObs = Predictability(dim=attractor$dim,theta=attractor$theta) +attractorMod<- CST_ProxiesAttractor(data=expL,quanti=qm) +predynMod = CST_Predictability(dim=attractorMod$dim,theta=attractorMod$theta) +dynbias= DynBiasCorrection(pred.dim.obs=predynObs$pred.dim,pred.theta.obs=predynMod$pred.theta, +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +scale atmospheric predictability.Nature Communications, 10(1), 1316. +DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +} diff --git a/man/Predictability.Rd b/man/Predictability.Rd new file mode 100644 index 00000000..2093a196 --- /dev/null +++ b/man/Predictability.Rd @@ -0,0 +1,72 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Predictability.R +\name{Predictability} +\alias{Predictability} +\title{Computing scores of predictability using two dinamical proxies +based on dynamical systems theory.} +\usage{ +Predictability(dim, theta) +} +\arguments{ +\item{dim}{data to create the attractor. Must be a matrix with the timesteps in nrow and the grids in ncol +for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints} + +\item{theta}{list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes} +} +\value{ +A list of length 2: +\itemize{ +\item\code{pred.dim} {a list of two lists 'qdim' and 'pos.d'. The 'qdim' list +contains values of local dimension 'dim' divided by terciles: +d1: lower tercile and more predictability, +d2: middle tercile, +d3: higher tercile and less predictability +The 'pos.d' list contains the position of each tercile in parameter 'dim'} + +\item\code{pred.theta} {a list of two lists 'qtheta' and 'pos.t'. +The 'qtheta' list contains values of the inverse of persistence 'theta' +divided by terciles: +th1: lower tercile and more predictability, +th2: middle tercile, +th3: higher tercile and less predictability +The 'pos.t' list contains the position of each tercile in parameter 'theta'} +} + +dyn_scores values from 0 to 1. A dyn_score of 1 indicates higher +predictability. +} +\description{ +This function divides in terciles the two dynamical proxies +computed with CST_ProxiesAttractor or ProxiesAttractor. These terciles will +be used to measure the predictability of the system in dyn_scores. When the +local dimension 'dim' is small and the inverse of persistence 'theta' is +small the predictability is higher, and viceversa. +} +\examples{ +# Creating an example of matrix dat(time,grids): +m <- matrix(rnorm(2000) * 10, nrow = 50, ncol = 40) +# imposing a threshold + quanti <- 0.90 +# computing dyn_scores from parameters dim and theta of the attractor +attractor <- ProxiesAttractor(dat = m, quanti = 0.60, iplot = FALSE) +predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +scale atmospheric predictability.Nature Communications, 10(1), 1316. +DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/ProxiesAttractor.Rd b/man/ProxiesAttractor.Rd new file mode 100644 index 00000000..497b9728 --- /dev/null +++ b/man/ProxiesAttractor.Rd @@ -0,0 +1,95 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_ProxiesAttractor.R, R/ProxiesAttractor.R +\name{ProxiesAttractor} +\alias{ProxiesAttractor} +\title{Computing two dinamical proxies of the attractor.} +\usage{ +ProxiesAttractor(data, quanti, iplot = FALSE, ncores = NULL) + +ProxiesAttractor(data, quanti, iplot = FALSE, ncores = NULL) +} +\arguments{ +\item{data}{data to create the attractor. Must be a matrix with the timesteps in nrow +and the grids in ncol(dat(time,grids)} + +\item{quanti}{list of arbitrary length of secondary grids. Each secondary grid is to +be provided as a list of length 2 with longitudes and latitudes} + +\item{iplot}{FALSE by default. If TRUE the function returns a simple plot +with three pannels on the top left a plot for local dimension 'd', on the top +right a plot for the inverse of the persistence 'theta', on the bottom a plot +for the 'attractor' plotted with the two properties 'd' (x axis) and 'theta' +(y axis)} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +dim and theta + +dim and theta +} +\description{ +This function computes two dinamical proxies of the attractor: +The local dimension (d) and the inverse of the persistence (theta). +These two parameters will be used as a condition for the computation of dynamical +scores to measure predictability and to compute bias correction conditioned by +the dynamics with the function DynBiasCorrection. +Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: + +This function computes two dinamical proxies of the attractor: +The local dimension (d) and the inverse of the persistence (theta). +These two parameters will be used as a condition for the computation of dynamical +scores to measure predictability and to compute bias correction conditioned by +the dynamics with the function DynBiasCorrection. +Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: +} +\examples{ +# Example 1: Computing the attractor using simple data +# Creating an example of matrix data(time,grids): +mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) +qm <- 0.90 # imposing a threshold +Attractor <- ProxiesAttractor(data = mat, quanti = qm) + +# Example 1: Computing the attractor using simple data +# Creating an example of matrix data(time,grids): +tm=2*6*3 # time +gm=4*10 # grids +m=matrix(rand(1,tm*gm),nrow=tm,ncol=gm) +qm=0.90 # imposing a threshold +Attractor=ProxiesAttractor(data=m,quanti=qm,iplot=FALSE) + +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. + +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., and Yiou, P. (2019). +The hammam effect or how a warm ocean enhances large scale atmospheric predictability. +Nature Communications, 10(1), 1316. DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} + +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} -- GitLab From c9260a6b2a35a61e27466ddd78af3bc8e8c74ef7 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Tue, 1 Jun 2021 13:14:16 +0000 Subject: [PATCH 10/19] Update DynBiasCorrection.R --- R/DynBiasCorrection.R | 184 +++++++++++++++++++++++++++++------------- 1 file changed, 127 insertions(+), 57 deletions(-) diff --git a/R/DynBiasCorrection.R b/R/DynBiasCorrection.R index f9bbd06c..8eda658c 100644 --- a/R/DynBiasCorrection.R +++ b/R/DynBiasCorrection.R @@ -3,16 +3,19 @@ #'properties of the data. #' #'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} #' #'@description This function perform a bias correction conditioned by the #'dynamical properties of the dataset. This function used the functions #''CST_Predictability' to divide in terciles the two dynamical proxies -#'computed with CST_ProxiesAttractor or ProxiesAttractor. A bias correction +#'computed with 'CST_ProxiesAttractor'. A bias correction #'between the model and the observations is performed using the division into #'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. #'For instance, model values with lower 'dim' will be corrected with observed #'values with lower 'dim', and the same for theta. The function gives two options -#'of bias correction: one for 'dim' and one for 'theta' +#'of bias correction: one for 'dim' and/or one for 'theta' #' #'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., #'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large @@ -22,69 +25,136 @@ #' Dynamical proxies of North Atlantic predictability and extremes. #' Scientific Reports, 7-41278, 2017. #' -#'@param pred.dim.mod output of CST_Predictability or a list of two lists 'qdim' and 'pos.d'. The 'qdim' list -#'contains values of local dimension 'dim' divided by terciles: -#'d1: lower tercile and more predictability, -#'d2: middle tercile, -#'d3: higher tercile and less predictability -#'The 'pos.d' list contains the position of each tercile in parameter 'dim' -#'@param pred.theta.mod output of CST_Predictability or a list of two lists 'qtheta' and 'pos.t'. -#'The 'qtheta' list contains values of the inverse of persistence 'theta' -#'divided by terciles: -#'th1: lower tercile and more predictability, -#'th2: middle tercile, -#'th3: higher tercile and less predictability -#'The 'pos.t' list contains the position of each tercile in parameter 'theta' -#' -#'@param pred.dim.obs same than pred.dim.obs but for the observations -#' -#'@param pred.theta.obs same than pred.theta.obs but for the observations -#' +#'@param obsL an s2vcube object with the observation data +#'@param expL an s2dvcube object with the experiment data +#'@param method the method to apply bias correction among these ones: +#'"PTF","RQUANT","QUANT","SSPLIN" +#'@param proxy the proxy (local dimension 'dim' or inverse of persistence +#''theta') to apply the dynamical conditioned bias correction method. +#'@param quanti quantile to perform the computation of local dimension and theta #'@param ncores The number of cores to use in parallel computation #' -#'@return dim.bias an s2dvcube object with a bias correction performed -#'conditioned by local dimension 'dim' -#' -#'@return theta.bias an s2dvcube object with a bias correction performed -#'conditioned by the inverse of the persistence 'theta' +#'@return dynbias an s2dvcube object with a bias correction performed +#'conditioned by local dimension 'dim' or inverse of persistence 'theta' #' +#'@import s2dverification +#'@import multiApply +#'@import qmap +#'@import CSTools #' #'@examples #'# example 1: simple data s2dvcube style -#' expL <- rnorm(1:200) -#' dim(expL) <- c(member=10,lat = 4, lon = 5) -#' obsL <- c(rnorm(1:180),expL[1,,]*1.2) -#' dim(obsL) <- c(time = 10,lat = 4, lon = 5) -#' time_obsL <- paste(rep("01", 50), rep("01", 50), 1953 : 2003, sep = "-") -#' time_expL <- time_obsL[1] -#' lon <- seq(-1,5,1.5) -#' lat <- seq(30,35,1.5) -#' qm=0.60 -#' expL <- s2dv_cube(data = expL, lat = lat, lon = lon, -#' Dates = list(start = time_expL, end = time_expL)) -#' obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, -#' Dates = list(start = time_obsL[1:10], end = time_obsL[1:10])) -#' attractor <- CST_ProxiesAttractor(data=obsL,quanti=qm) -#' predynObs = Predictability(dim=attractor$dim,theta=attractor$theta) -#' attractorMod<- CST_ProxiesAttractor(data=expL,quanti=qm) -#' predynMod = CST_Predictability(dim=attractorMod$dim,theta=attractorMod$theta) -#' dynbias= DynBiasCorrection(pred.dim.obs=predynObs$pred.dim,pred.theta.obs=predynMod$pred.theta, -# pred.dim.mod=predynObs$pred.dim,pred.theta.mod=predynMod$pred.theta) -#'@export -DynBiasCorrection<- function(pred.dim.obs,pred.theta.obs,pred.dim.mod, - pred.theta.mod,ncores=NULL){ - if (is.null(pred.dim.mod)) { - stop("Parameter 'pred.dim.mod' is mandatory") +#'expL <- rnorm(1:2000) +#'dim (expL) <- c(time =100,lat = 4, lon = 5) +#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +#'dim (obsL) <- c(time = 100,lat = 4, lon = 5) +#'time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") +#'time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") +#'lon <- seq(-1,5,1.5) +#'lat <- seq(30,35,1.5) +#'# qm=0.98 # too high for this short dataset, it is possible that doesn't +#'# get the requirement, in that case it would be necessary select a lower qm +#'# for instance qm=0.60 +#'qm=0.60 +#'method="QUANT" +#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' Dates = list(start = time_expL, end = time_expL)) +#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' Dates = list(start = time_obsL, end = time_obsL)) +#' +#'dynbias<- DynBiasCorrection(obsL = obsL, expL = expL, method = method, +#' proxy= "dim", quanti=qm, ncores=NULL) +#' + +DynBiasCorrection<- function(obsL = obsL, expL = expL, + method = method, + proxy= "dim", quanti=qm, ncores=NULL){ + if (is.null(obsL)) { + stop("Parameter 'obsL' is mandatory") } - if (is.null(pred.dim.obs)) { - stop("Parameter 'pred.dim.obs' is mandatory") + if (is.null(expL)) { + stop("Parameter 'expL' is mandatory") } - if (is.null(pred.theta.mod)) { - stop("Parameter 'pred.dim.mod' is mandatory") + if (is.null(method)) { + stop("Parameter 'method' is mandatory") } - if (is.null(pred.theta.obs)) { - stop("Parameter 'pred.dim.obs' is mandatory") + if (is.null(quanti)) { + stop("Parameter 'quanti' is mandatory") + } + if (is.null(proxy)) { + stop("Parameter 'proxy' is mandatory") + } + if (!inherits(obsL, 's2dv_cube')) { + stop("Parameter 'obsL' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!inherits(expL, 's2dv_cube')) { + stop("Parameter 'expL' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + + attractor.obs <- CST_ProxiesAttractor(data=obsL,quanti=qm) + predyn.obs <- CST_Predictability(dim=attractor.obs$dim, + theta=attractor.obs$theta) + attractor.exp <- CST_ProxiesAttractor(data=expL,quanti=qm) + predyn.exp <- CST_Predictability(dim=attractor.exp$dim, + theta=attractor.exp$theta) + + fun1 <- function(data,list1,tar_dim){ + lapply(list1, function(x){ + Subset(data,indices=x,along=tar_dim)})} + + if(proxy=="dim"){ + obs= fun1(data=obsL$data,list1=predyn.obs$pred.dim$pos.d,tar_dim='time') + exp= fun1(data=expL$data,list1=predyn.exp$pred.dim$pos.d,tar_dim='time') + } + if(proxy=="theta"){ + obs= fun1(data=obsL$data,list1=predyn.obs$pred.theta$pos.t,tar_dim='time') + exp= fun1(data=expL$data,list1=predyn.exp$pred.theta$pos.t,tar_dim='time') + } + + biascorrection <- function(obs,exp,method){ + ## functions fitQmap and doQmap + if(method=="PTF"){ + qm.fit <- fitQmap(obs,exp, + method="PTF", + transfun="expasympt", + cost="RSS",wett.day=TRUE) + qmap <- doQmap(exp,qm.fit)} + if(method=="QUANT"){ + qm.fit <- fitQmap(obs,exp, + method="QUANT",qstep=0.01) + qmap <- doQmap(exp,qm.fit,type="tricub")} + if(method=="RQUANT"){ + qm.fit <- fitQmap(obs,exp, + method="RQUANT",qstep=0.01) + qmap <- doQmap(exp,qm.fit,type="linear")} + if(method=="SSPLIN"){ + qm.fit <- fitQmap(obs,exp,qstep=0.01, + method="SSPLIN") + qmap <- doQmap(exp,qm.fit)} + + return(qmap) + } + ## functions bias correction CSTools + # ? + latlonApply <- function(obs,exp,method){ + Apply(list(exp=exp,obs=obs),target_dims='time', + fun=biascorrection,method=method)} + + qmap <- mapply(latlonApply,obs,exp,method) + expAdjust = array(dim=dim(expL$data)) + + if(proxy=="dim"){ + expAdjust[predyn.exp$pred.dim$pos.d$d1,,]=qmap$d1.output1 + expAdjust[predyn.exp$pred.dim$pos.d$d2,,]=qmap$d2.output1 + expAdjust[predyn.exp$pred.dim$pos.d$d3,,]=qmap$d3.output1 + } + if(proxy=="theta"){ + expAdjust[predyn.exp$pred.theta$pos.t$th1,,]=qmap$th1.output1 + expAdjust[predyn.exp$pred.theta$pos.t$th2,,]=qmap$th2.output1 + expAdjust[predyn.exp$pred.theta$pos.t$th3,,]=qmap$th3.output1 } - dim.bias <- theta.bis <- NULL - return(list(dim.bias=dim.bias,theta.bias=theta.bias)) + + return(list(expAdjust=expAdjust)) } -- GitLab From 3c7e14914f343b2a100758689e2dbd1750a47bf2 Mon Sep 17 00:00:00 2001 From: nperez Date: Wed, 2 Jun 2021 11:42:02 +0200 Subject: [PATCH 11/19] DynBiasCorrection using atomic function --- NAMESPACE | 1 + R/CST_DynBiasCorrection.R | 202 +++++++++++++++++++++++++++++++++++ R/DynBiasCorrection.R | 160 --------------------------- R/Predictability.R | 1 + man/CST_DynBiasCorrection.Rd | 88 +++++++++++++++ man/DynBiasCorrection.Rd | 77 ++++++------- man/Predictability.Rd | 1 + 7 files changed, 324 insertions(+), 206 deletions(-) create mode 100644 R/CST_DynBiasCorrection.R delete mode 100644 R/DynBiasCorrection.R create mode 100644 man/CST_DynBiasCorrection.Rd diff --git a/NAMESPACE b/NAMESPACE index f2da859a..26042662 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,6 +12,7 @@ export(CST_BEI_Weighting) export(CST_BiasCorrection) export(CST_Calibration) export(CST_CategoricalEnsCombination) +export(CST_DynBiasCorrection) export(CST_EnsClustering) export(CST_Load) export(CST_MergeDims) diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R new file mode 100644 index 00000000..7149e49f --- /dev/null +++ b/R/CST_DynBiasCorrection.R @@ -0,0 +1,202 @@ +#'@rdname CST_DynBiasCorrection +#'@title Performing a Bias Correction conditioned by the dynamical +#'properties of the data. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function perform a bias correction conditioned by the +#'dynamical properties of the dataset. This function internally uses the functions +#''Predictability' to divide in terciles the two dynamical proxies +#'computed with 'CST_ProxiesAttractor'. A bias correction +#'between the model and the observations is performed using the division into +#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +#'For instance, model values with lower 'dim' will be corrected with observed +#'values with lower 'dim', and the same for theta. The function gives two options +#'of bias correction: one for 'dim' and/or one for 'theta' +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param exp an s2v_cube object with the experiment data +#'@param obs an s2dv_cube object with the reference data +#'@param method a character string indicating the method to apply bias correction among these ones: +#'"PTF","RQUANT","QUANT","SSPLIN" +#'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. +#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta +#'@param time_dim a character string indicating the name of the temporal dimension +#'@param ncores The number of cores to use in parallel computation +#' +#'@return dynbias an s2dvcube object with a bias correction performed +#'conditioned by local dimension 'dim' or inverse of persistence 'theta' +#' +#'@examples +#'# example 1: simple data s2dvcube style +#'set.seed(1) +#'expL <- rnorm(1:2000) +#'dim (expL) <- c(time =100,lat = 4, lon = 5) +#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +#'dim (obsL) <- c(time = 100,lat = 4, lon = 5) +#'time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") +#'time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") +#'lon <- seq(-1,5,1.5) +#'lat <- seq(30,35,1.5) +# qm=0.98 # too high for this short dataset, it is possible that doesn't +#'# get the requirement, in that case it would be necessary select a lower qm +#'# for instance qm=0.60 +#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' Dates = list(start = time_expL, end = time_expL)) +#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' Dates = list(start = time_obsL, end = time_obsL)) +#'dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", +#' quanti = 0.6, time_dim = 'time') +#' +#'@export +CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', + proxy = "dim", quanti, time_dim = 'ftime', + ncores = NULL) { + if (!inherits(obs, 's2dv_cube')) { + stop("Parameter 'obs' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + if (!inherits(exp, 's2dv_cube')) { + stop("Parameter 'exp' must be of the class 's2dv_cube', ", + "as output by CSTools::CST_Load.") + } + exp$data <- DynBiasCorrection(exp = exp$data, obs = obs$data, method = method, + proxy = proxy, quanti = quanti, + time_dim = time_dim, ncores = ncores) + return(exp) +} +#'@rdname DynBiasCorrection +#'@title Performing a Bias Correction conditioned by the dynamical +#'properties of the data. +#' +#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} +#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} +#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} +#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +#' +#'@description This function perform a bias correction conditioned by the +#'dynamical properties of the dataset. This function used the functions +#''CST_Predictability' to divide in terciles the two dynamical proxies +#'computed with 'CST_ProxiesAttractor'. A bias correction +#'between the model and the observations is performed using the division into +#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +#'For instance, model values with lower 'dim' will be corrected with observed +#'values with lower 'dim', and the same for theta. The function gives two options +#'of bias correction: one for 'dim' and/or one for 'theta' +#' +#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +#'scale atmospheric predictability.Nature Communications, 10(1), 1316. +#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " +#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +#' Dynamical proxies of North Atlantic predictability and extremes. +#' Scientific Reports, 7-41278, 2017. +#' +#'@param exp a multidimensional array with named dimensions with the experiment data +#'@param obs a multidimensional array with named dimensions with the observation data +#'@param method a character string indicating the method to apply bias correction among these ones: +#'"PTF","RQUANT","QUANT","SSPLIN" +#'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. +#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta +#'@param time_dim a character string indicating the name of the temporal dimension +#'@param ncores The number of cores to use in parallel computation +#' +#'@return a multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta' +#' +#'@import multiApply +#'@importFrom s2dverification Subset +#'@import qmap +#'@examples +#'expL <- rnorm(1:2000) +#'dim (expL) <- c(time =100,lat = 4, lon = 5) +#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +#'dim (obsL) <- c(time = 100,lat = 4, lon = 5) +#'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, +#' proxy= "dim", quanti = 0.6) +#'@export +DynBiasCorrection<- function(exp, obs, method = 'QUANT', + proxy = "dim", quanti, + time_dim = 'time', ncores = NULL){ + if (is.null(obs)) { + stop("Parameter 'obs' cannot be NULL.") + } + if (is.null(exp)) { + stop("Parameter 'exp' cannot be NULL.") + } + if (is.null(method)) { + stop("Parameter 'method' cannot be NULL.") + } + if (is.null(quanti)) { + stop("Parameter 'quanti' cannot be NULL.") + } + if (is.null(proxy)) { + stop("Parameter 'proxy' cannot be NULL.") + } + dims <- dim(exp) + + attractor.obs <- ProxiesAttractor(data = obs, quanti = quanti) + predyn.obs <- Predictability(dim = attractor.obs$dim, + theta = attractor.obs$theta) + attractor.exp <- ProxiesAttractor(data = exp, quanti = quanti) + predyn.exp <- Predictability(dim = attractor.exp$dim, + theta = attractor.exp$theta) + + if (proxy == "dim") { + adjusted <- Apply(list(exp, obs), target_dims = time_dim, + fun = .dynbias, method, + predyn.exp = predyn.exp$pred.dim$pos.d, + predyn.obs = predyn.obs$pred.dim$pos.d, + ncores = ncores, output_dims = time_dim)$output1 + } else if (proxy == "theta") { + adjusted <- Apply(list(exp, obs), target_dims = time_dim, + fun = .dynbias, method, + predyn.exp = predyn.exp$pred.theta$pos.t, + predyn.obs = predyn.obs$pred.theta$pos.t, + ncores = ncores, output_dims = time_dim)$output1 + } else { + stop ("Parameter 'proxy' must be set as 'dim' or 'theta'.") + } + return(adjusted) +} + +.dynbias <- function(exp, obs, method, predyn.exp, predyn.obs) { + result <- array(rep(NA, length(exp))) + res <- lapply(1:3, function(x) { + exp_sub <- exp[predyn.exp[[x]]] + obs_sub <- obs[predyn.obs[[x]]] + adjust <- .qbiascorrection(exp_sub, obs_sub, method) + result[predyn.exp[[x]]] <<- adjust + return(NULL) + }) + return(result) +} +.qbiascorrection <- function(expX, obsX, method) { + ## functions fitQmap and doQmap + if (method == "PTF") { + qm.fit <- fitQmap(obsX, expX, method = "PTF", transfun = "expasympt", + cost = "RSS", wett.day = TRUE) + qmap <- doQmap(expX, qm.fit) + } else if (method == "QUANT") { + qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01) + qmap <- doQmap(expX, qm.fit, type = "tricub") + } else if (method == "RQUANT") { + qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01) + qmap <- doQmap(expX, qm.fit, type = "linear") + } else if (method == "SSPLIN") { + qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN") + qmap <- doQmap(expX, qm.fit) + } else { + stop ("Parameter 'method' doesn't match any of the available methods.") + } + return(qmap) +} diff --git a/R/DynBiasCorrection.R b/R/DynBiasCorrection.R deleted file mode 100644 index 8eda658c..00000000 --- a/R/DynBiasCorrection.R +++ /dev/null @@ -1,160 +0,0 @@ -#'@rdname DynBiasCorrection -#'@title Performing a Bias Correction conditioned by the dynamical -#'properties of the data. -#' -#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} -#'@author Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} -#'@author Veronica Torralba, \email{veronica.torralba@cmcc.it} -#'@author Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} -#' -#'@description This function perform a bias correction conditioned by the -#'dynamical properties of the dataset. This function used the functions -#''CST_Predictability' to divide in terciles the two dynamical proxies -#'computed with 'CST_ProxiesAttractor'. A bias correction -#'between the model and the observations is performed using the division into -#'terciles of the local dimension 'dim' and inverse of the persistence 'theta'. -#'For instance, model values with lower 'dim' will be corrected with observed -#'values with lower 'dim', and the same for theta. The function gives two options -#'of bias correction: one for 'dim' and/or one for 'theta' -#' -#'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., -#'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large -#'scale atmospheric predictability.Nature Communications, 10(1), 1316. -#'DOI = https://doi.org/10.1038/s41467-019-09305-8 " -#'@references Faranda, D., Gabriele Messori and Pascal Yiou. (2017). -#' Dynamical proxies of North Atlantic predictability and extremes. -#' Scientific Reports, 7-41278, 2017. -#' -#'@param obsL an s2vcube object with the observation data -#'@param expL an s2dvcube object with the experiment data -#'@param method the method to apply bias correction among these ones: -#'"PTF","RQUANT","QUANT","SSPLIN" -#'@param proxy the proxy (local dimension 'dim' or inverse of persistence -#''theta') to apply the dynamical conditioned bias correction method. -#'@param quanti quantile to perform the computation of local dimension and theta -#'@param ncores The number of cores to use in parallel computation -#' -#'@return dynbias an s2dvcube object with a bias correction performed -#'conditioned by local dimension 'dim' or inverse of persistence 'theta' -#' -#'@import s2dverification -#'@import multiApply -#'@import qmap -#'@import CSTools -#' -#'@examples -#'# example 1: simple data s2dvcube style -#'expL <- rnorm(1:2000) -#'dim (expL) <- c(time =100,lat = 4, lon = 5) -#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) -#'dim (obsL) <- c(time = 100,lat = 4, lon = 5) -#'time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") -#'time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") -#'lon <- seq(-1,5,1.5) -#'lat <- seq(30,35,1.5) -#'# qm=0.98 # too high for this short dataset, it is possible that doesn't -#'# get the requirement, in that case it would be necessary select a lower qm -#'# for instance qm=0.60 -#'qm=0.60 -#'method="QUANT" -#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, -#' Dates = list(start = time_expL, end = time_expL)) -#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, -#' Dates = list(start = time_obsL, end = time_obsL)) -#' -#'dynbias<- DynBiasCorrection(obsL = obsL, expL = expL, method = method, -#' proxy= "dim", quanti=qm, ncores=NULL) -#' - -DynBiasCorrection<- function(obsL = obsL, expL = expL, - method = method, - proxy= "dim", quanti=qm, ncores=NULL){ - if (is.null(obsL)) { - stop("Parameter 'obsL' is mandatory") - } - if (is.null(expL)) { - stop("Parameter 'expL' is mandatory") - } - if (is.null(method)) { - stop("Parameter 'method' is mandatory") - } - if (is.null(quanti)) { - stop("Parameter 'quanti' is mandatory") - } - if (is.null(proxy)) { - stop("Parameter 'proxy' is mandatory") - } - if (!inherits(obsL, 's2dv_cube')) { - stop("Parameter 'obsL' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.") - } - if (!inherits(expL, 's2dv_cube')) { - stop("Parameter 'expL' must be of the class 's2dv_cube', ", - "as output by CSTools::CST_Load.") - } - - attractor.obs <- CST_ProxiesAttractor(data=obsL,quanti=qm) - predyn.obs <- CST_Predictability(dim=attractor.obs$dim, - theta=attractor.obs$theta) - attractor.exp <- CST_ProxiesAttractor(data=expL,quanti=qm) - predyn.exp <- CST_Predictability(dim=attractor.exp$dim, - theta=attractor.exp$theta) - - fun1 <- function(data,list1,tar_dim){ - lapply(list1, function(x){ - Subset(data,indices=x,along=tar_dim)})} - - if(proxy=="dim"){ - obs= fun1(data=obsL$data,list1=predyn.obs$pred.dim$pos.d,tar_dim='time') - exp= fun1(data=expL$data,list1=predyn.exp$pred.dim$pos.d,tar_dim='time') - } - if(proxy=="theta"){ - obs= fun1(data=obsL$data,list1=predyn.obs$pred.theta$pos.t,tar_dim='time') - exp= fun1(data=expL$data,list1=predyn.exp$pred.theta$pos.t,tar_dim='time') - } - - biascorrection <- function(obs,exp,method){ - ## functions fitQmap and doQmap - if(method=="PTF"){ - qm.fit <- fitQmap(obs,exp, - method="PTF", - transfun="expasympt", - cost="RSS",wett.day=TRUE) - qmap <- doQmap(exp,qm.fit)} - if(method=="QUANT"){ - qm.fit <- fitQmap(obs,exp, - method="QUANT",qstep=0.01) - qmap <- doQmap(exp,qm.fit,type="tricub")} - if(method=="RQUANT"){ - qm.fit <- fitQmap(obs,exp, - method="RQUANT",qstep=0.01) - qmap <- doQmap(exp,qm.fit,type="linear")} - if(method=="SSPLIN"){ - qm.fit <- fitQmap(obs,exp,qstep=0.01, - method="SSPLIN") - qmap <- doQmap(exp,qm.fit)} - - return(qmap) - } - ## functions bias correction CSTools - # ? - latlonApply <- function(obs,exp,method){ - Apply(list(exp=exp,obs=obs),target_dims='time', - fun=biascorrection,method=method)} - - qmap <- mapply(latlonApply,obs,exp,method) - expAdjust = array(dim=dim(expL$data)) - - if(proxy=="dim"){ - expAdjust[predyn.exp$pred.dim$pos.d$d1,,]=qmap$d1.output1 - expAdjust[predyn.exp$pred.dim$pos.d$d2,,]=qmap$d2.output1 - expAdjust[predyn.exp$pred.dim$pos.d$d3,,]=qmap$d3.output1 - } - if(proxy=="theta"){ - expAdjust[predyn.exp$pred.theta$pos.t$th1,,]=qmap$th1.output1 - expAdjust[predyn.exp$pred.theta$pos.t$th2,,]=qmap$th2.output1 - expAdjust[predyn.exp$pred.theta$pos.t$th3,,]=qmap$th3.output1 - } - - return(list(expAdjust=expAdjust)) -} diff --git a/R/Predictability.R b/R/Predictability.R index 1989ef3d..701b812e 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -49,6 +49,7 @@ #'@examples #'# Creating an example of matrix dat(time,grids): #'m <- matrix(rnorm(2000) * 10, nrow = 50, ncol = 40) +#'names(dim(m)) <- c('time', 'grid') #'# imposing a threshold #' quanti <- 0.90 #'# computing dyn_scores from parameters dim and theta of the attractor diff --git a/man/CST_DynBiasCorrection.Rd b/man/CST_DynBiasCorrection.Rd new file mode 100644 index 00000000..facf6f51 --- /dev/null +++ b/man/CST_DynBiasCorrection.Rd @@ -0,0 +1,88 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/CST_DynBiasCorrection.R +\name{CST_DynBiasCorrection} +\alias{CST_DynBiasCorrection} +\title{Performing a Bias Correction conditioned by the dynamical +properties of the data.} +\usage{ +CST_DynBiasCorrection( + exp, + obs, + method = "QUANT", + proxy = "dim", + quanti, + time_dim = "ftime", + ncores = NULL +) +} +\arguments{ +\item{exp}{an s2v_cube object with the experiment data} + +\item{obs}{an s2dv_cube object with the reference data} + +\item{method}{a character string indicating the method to apply bias correction among these ones: +"PTF","RQUANT","QUANT","SSPLIN"} + +\item{proxy}{a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} + +\item{time_dim}{a character string indicating the name of the temporal dimension} + +\item{ncores}{The number of cores to use in parallel computation} +} +\value{ +dynbias an s2dvcube object with a bias correction performed +conditioned by local dimension 'dim' or inverse of persistence 'theta' +} +\description{ +This function perform a bias correction conditioned by the +dynamical properties of the dataset. This function internally uses the functions +'Predictability' to divide in terciles the two dynamical proxies +computed with 'CST_ProxiesAttractor'. A bias correction +between the model and the observations is performed using the division into +terciles of the local dimension 'dim' and inverse of the persistence 'theta'. +For instance, model values with lower 'dim' will be corrected with observed +values with lower 'dim', and the same for theta. The function gives two options +of bias correction: one for 'dim' and/or one for 'theta' +} +\examples{ +# example 1: simple data s2dvcube style +set.seed(1) +expL <- rnorm(1:2000) +dim (expL) <- c(time =100,lat = 4, lon = 5) +obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +dim (obsL) <- c(time = 100,lat = 4, lon = 5) +time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") +time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") +lon <- seq(-1,5,1.5) +lat <- seq(30,35,1.5) +# get the requirement, in that case it would be necessary select a lower qm +# for instance qm=0.60 +expL <- s2dv_cube(data = expL, lat = lat, lon = lon, + Dates = list(start = time_expL, end = time_expL)) +obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, + Dates = list(start = time_obsL, end = time_obsL)) +dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", + quanti = 0.6, time_dim = 'time') + +} +\references{ +Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., +and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large +scale atmospheric predictability.Nature Communications, 10(1), 1316. +DOI = https://doi.org/10.1038/s41467-019-09305-8 " + +Faranda, D., Gabriele Messori and Pascal Yiou. (2017). +Dynamical proxies of North Atlantic predictability and extremes. +Scientific Reports, 7-41278, 2017. +} +\author{ +Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} +} diff --git a/man/DynBiasCorrection.Rd b/man/DynBiasCorrection.Rd index 878fc764..bd60e0f3 100644 --- a/man/DynBiasCorrection.Rd +++ b/man/DynBiasCorrection.Rd @@ -1,78 +1,57 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/DynBiasCorrection.R +% Please edit documentation in R/CST_DynBiasCorrection.R \name{DynBiasCorrection} \alias{DynBiasCorrection} \title{Performing a Bias Correction conditioned by the dynamical properties of the data.} \usage{ DynBiasCorrection( - pred.dim.obs, - pred.theta.obs, - pred.dim.mod, - pred.theta.mod, + exp, + obs, + method = "QUANT", + proxy = "dim", + quanti, + time_dim = "time", ncores = NULL ) } \arguments{ -\item{pred.dim.obs}{same than pred.dim.obs but for the observations} +\item{exp}{a multidimensional array with named dimensions with the experiment data} -\item{pred.theta.obs}{same than pred.theta.obs but for the observations} +\item{obs}{a multidimensional array with named dimensions with the observation data} -\item{pred.dim.mod}{output of CST_Predictability or a list of two lists 'qdim' and 'pos.d'. The 'qdim' list -contains values of local dimension 'dim' divided by terciles: -d1: lower tercile and more predictability, -d2: middle tercile, -d3: higher tercile and less predictability -The 'pos.d' list contains the position of each tercile in parameter 'dim'} +\item{method}{a character string indicating the method to apply bias correction among these ones: +"PTF","RQUANT","QUANT","SSPLIN"} -\item{pred.theta.mod}{output of CST_Predictability or a list of two lists 'qtheta' and 'pos.t'. -The 'qtheta' list contains values of the inverse of persistence 'theta' -divided by terciles: -th1: lower tercile and more predictability, -th2: middle tercile, -th3: higher tercile and less predictability -The 'pos.t' list contains the position of each tercile in parameter 'theta'} +\item{proxy}{a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} + +\item{time_dim}{a character string indicating the name of the temporal dimension} \item{ncores}{The number of cores to use in parallel computation} } \value{ -dim.bias an s2dvcube object with a bias correction performed -conditioned by local dimension 'dim' - -theta.bias an s2dvcube object with a bias correction performed -conditioned by the inverse of the persistence 'theta' +a multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta' } \description{ This function perform a bias correction conditioned by the dynamical properties of the dataset. This function used the functions 'CST_Predictability' to divide in terciles the two dynamical proxies -computed with CST_ProxiesAttractor or ProxiesAttractor. A bias correction +computed with 'CST_ProxiesAttractor'. A bias correction between the model and the observations is performed using the division into terciles of the local dimension 'dim' and inverse of the persistence 'theta'. For instance, model values with lower 'dim' will be corrected with observed values with lower 'dim', and the same for theta. The function gives two options -of bias correction: one for 'dim' and one for 'theta' +of bias correction: one for 'dim' and/or one for 'theta' } \examples{ -# example 1: simple data s2dvcube style -expL <- rnorm(1:200) -dim(expL) <- c(member=10,lat = 4, lon = 5) -obsL <- c(rnorm(1:180),expL[1,,]*1.2) -dim(obsL) <- c(time = 10,lat = 4, lon = 5) -time_obsL <- paste(rep("01", 50), rep("01", 50), 1953 : 2003, sep = "-") -time_expL <- time_obsL[1] -lon <- seq(-1,5,1.5) -lat <- seq(30,35,1.5) -qm=0.60 -expL <- s2dv_cube(data = expL, lat = lat, lon = lon, - Dates = list(start = time_expL, end = time_expL)) -obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, - Dates = list(start = time_obsL[1:10], end = time_obsL[1:10])) -attractor <- CST_ProxiesAttractor(data=obsL,quanti=qm) -predynObs = Predictability(dim=attractor$dim,theta=attractor$theta) -attractorMod<- CST_ProxiesAttractor(data=expL,quanti=qm) -predynMod = CST_Predictability(dim=attractorMod$dim,theta=attractorMod$theta) -dynbias= DynBiasCorrection(pred.dim.obs=predynObs$pred.dim,pred.theta.obs=predynMod$pred.theta, +expL <- rnorm(1:2000) +dim (expL) <- c(time =100,lat = 4, lon = 5) +obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +dim (obsL) <- c(time = 100,lat = 4, lon = 5) +dynbias <- DynBiasCorrection(exp = expL, obs = obsL, + proxy= "dim", quanti = 0.6) } \references{ Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., @@ -86,4 +65,10 @@ Scientific Reports, 7-41278, 2017. } \author{ Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} + +Maria M. Chaves-Montero, \email{mdm.chaves-montero@cmcc.it} + +Veronica Torralba, \email{veronica.torralba@cmcc.it} + +Davide Faranda, \email{davide.faranda@lsce.ipsl.fr} } diff --git a/man/Predictability.Rd b/man/Predictability.Rd index 2093a196..be0fb623 100644 --- a/man/Predictability.Rd +++ b/man/Predictability.Rd @@ -45,6 +45,7 @@ small the predictability is higher, and viceversa. \examples{ # Creating an example of matrix dat(time,grids): m <- matrix(rnorm(2000) * 10, nrow = 50, ncol = 40) +names(dim(m)) <- c('time', 'grid') # imposing a threshold quanti <- 0.90 # computing dyn_scores from parameters dim and theta of the attractor -- GitLab From 8cb1ef9b8c982e63cf742a109e9cb49ece043e6c Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Mon, 7 Jun 2021 19:46:35 +0000 Subject: [PATCH 12/19] corrections review CST_ProxiesAttractor.R --- R/CST_ProxiesAttractor.R | 42 ++++++++++++++-------------------------- 1 file changed, 15 insertions(+), 27 deletions(-) diff --git a/R/CST_ProxiesAttractor.R b/R/CST_ProxiesAttractor.R index 7df33fee..518968fa 100644 --- a/R/CST_ProxiesAttractor.R +++ b/R/CST_ProxiesAttractor.R @@ -23,8 +23,7 @@ #'@param data a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow #'and the grids in ncol(dat(time,grids) # -#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to -#' be provided as a list of length 2 with longitudes and latitudes +#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta #' #'@param ncores The number of cores to use in parallel computation #' @@ -70,18 +69,11 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ #' Scientific Reports, 7-41278, 2017. #' #'@param data a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'. -#'@param quanti list of arbitrary length of secondary grids. Each secondary grid is to -#' be provided as a list of length 2 with longitudes and latitudes. -#'@param iplot a logical indicating if to plot results. +#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta #'@param ncores The number of cores to use in parallel computation. -#' -#'@param iplot FALSE by default. If TRUE the function returns a simple plot -#'with three pannels on the top left a plot for local dimension 'd', on the top -#'right a plot for the inverse of the persistence 'theta', on the bottom a plot -#'for the 'attractor' plotted with the two properties 'd' (x axis) and 'theta' -#'(y axis) #' #'@return dim and theta +#' #'@import multiApply #' #'@examples @@ -90,18 +82,24 @@ CST_ProxiesAttractor <- function(data, quanti, ncores = NULL){ #'mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) #'qm <- 0.90 # imposing a threshold #'Attractor <- ProxiesAttractor(data = mat, quanti = qm) +#'# to plot the result +#'time = c(1:length(Attractor$theta)) +#'layout(matrix(c(1, 3, 2, 3), 2, 2)) +#'plot(time, Attractor$dim, xlab = 'time', ylab = 'd', +#' main = 'local dimension', type = 'l') +#'plot(time, Attractor$theta, xlab = 'time', ylab = 'theta', main = 'theta') +#'plot(Attractor$dim, Attractor$theta, col = 'blue', +#' main = "Proxies of the Attractor", +#' xlab = "local dimension", ylab = "theta", lwd = 8, 'p') #' #'@export -ProxiesAttractor <- function(data, quanti, iplot = FALSE, ncores = NULL){ +ProxiesAttractor <- function(data, quanti, ncores = NULL){ if (is.null(data)) { stop("Parameter 'data' cannot be NULL.") } if (is.null(quanti)) { stop("Parameter 'quanti' is mandatory") - } - if (!is.logical(iplot) || is.null(iplot)) { - stop("Parameter 'iplot' is required and needs to be TRUE or FALSE.") } if (any(names(dim(data)) %in% 'sdate')) { if (any(names(dim(data)) %in% 'ftime')) { @@ -129,7 +127,7 @@ ProxiesAttractor <- function(data, quanti, iplot = FALSE, ncores = NULL){ } attractor <- Apply(data, target_dims = c('time', 'grid'), fun = .proxiesattractor, - quanti = quanti , iplot = FALSE, ncores = ncores) + quanti = quanti , ncores = ncores) # rename dimensions attractor <- lapply(attractor, FUN = function(x, dimname){ @@ -139,7 +137,7 @@ ProxiesAttractor <- function(data, quanti, iplot = FALSE, ncores = NULL){ return(list(dim = attractor$dim, theta = attractor$theta)) } -.proxiesattractor <- function(data, quanti, iplot = FALSE) { +.proxiesattractor <- function(data, quanti) { # expected dimensions data: time and grid logdista <- Apply(data, target_dims = 'grid', fun = function(x, y){ @@ -176,16 +174,6 @@ ProxiesAttractor <- function(data, quanti, iplot = FALSE, ncores = NULL){ names(dim(logdista)) <- c('dim1', 'dim2') proxies <- Apply(data = list(logdista = logdista), target_dims = list('dim1'), fun = Theta, quanti = quanti) - if(iplot == TRUE) { - time = c(1:length(proxies$theta)) - layout(matrix(c(1, 3, 2, 3), 2, 2)) - plot(time, proxies$dim, xlab = 'time', ylab = 'd', - main = 'local dimension', type = 'l') - plot(time, proxies$theta, xlab = 'time', ylab = 'theta', main = 'theta') - plot(proxies$dim, proxies$theta, col = 'blue', - main = "Proxies of the Attractor", - xlab = "local dimension", ylab = "theta", lwd = 8, 'p') - } return(list(dim = proxies$dim, theta = proxies$theta)) } -- GitLab From bdd294a0a68f5ede446f620afd8658af7e869e95 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Mon, 7 Jun 2021 19:48:17 +0000 Subject: [PATCH 13/19] Corrections review Predictability.R --- R/Predictability.R | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/R/Predictability.R b/R/Predictability.R index 701b812e..d17ea24c 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -56,7 +56,28 @@ #' attractor <- ProxiesAttractor(dat = m, quanti = 0.60, iplot = FALSE) #' predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) #'@export -Predictability <- function(dim, theta) { +#' +Predictability<- function(dim, theta, ncores = NULL) { + # if (!inherits(dim, 's2dv_cube')) { + # stop("Parameter 'dim' must be of the class 's2dv_cube', ", + # "as output by CSTools::CST_Load.") + # } + # if (!inherits(theta, 's2dv_cube')) { + # stop("Parameter 'theta' must be of the class 's2dv_cube', ", + # "as output by CSTools::CST_Load.") + # } + + pred <- Apply(list(dim, theta), target_dims = 'time', + fun = .predictability, + ncores = ncores) + dim(pred$dyn_scores) <- dim(theta) + return(list(pred.dim = list(qdim = list(pred$qdim.d1,pred$qdim.d2,pred$qdim.d3), + pos.d = list(pred$pos.d1,pred$pos.d2,pred$pos.d3)), + pred.theta = list(qtheta = list(pred$qtheta.th1,pred$qtheta.th2,pred$qtheta.th3), + pos.t = list(pred$pos.th1,pred$pos.th2,pred$pos.th3)), + dyn_scores = pred$dyn_scores)) +} +.predictability <- function(dim, theta) { if (is.null(dim)) { stop("Parameter 'dim' cannot be NULL.") } @@ -108,9 +129,13 @@ Predictability <- function(dim, theta) { dyn_scores[which(pos %in% d3th2)]<- 2/9 dyn_scores[which(pos %in% d3th3)]<- 1/9 -return(list(pred.dimi =list(qdim = qdim, pos.d = pos.d), - pred.theta = list(qtheta = qtheta, pos.t = pos.t), +return(list(qdim.d1 = dim[d1], qdim.d2 = dim[d2], qdim.d3 = dim[d3], + pos.d1 = d1, pos.d2 = d2, pos.d3 =d3, + qtheta.th1 = theta[th1], qtheta.th2 = theta[th2], qtheta.th3 = theta[th3], + pos.th1 = th1, pos.th2 = th2, pos.th3 = th3, dyn_scores = dyn_scores)) } + + -- GitLab From 4efac0fbb14ff551f4800e2c1a615cde94468120 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Mon, 7 Jun 2021 19:49:50 +0000 Subject: [PATCH 14/19] Corrections review in CST_DynBiasCorrection.R --- R/CST_DynBiasCorrection.R | 61 ++++++++++++++++++++++++++++++++------- 1 file changed, 50 insertions(+), 11 deletions(-) diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R index 7149e49f..00d8dc68 100644 --- a/R/CST_DynBiasCorrection.R +++ b/R/CST_DynBiasCorrection.R @@ -60,7 +60,7 @@ #' #'@export CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', - proxy = "dim", quanti, time_dim = 'ftime', + proxy = "dim", quanti, ncores = NULL) { if (!inherits(obs, 's2dv_cube')) { stop("Parameter 'obs' must be of the class 's2dv_cube', ", @@ -71,8 +71,7 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', "as output by CSTools::CST_Load.") } exp$data <- DynBiasCorrection(exp = exp$data, obs = obs$data, method = method, - proxy = proxy, quanti = quanti, - time_dim = time_dim, ncores = ncores) + proxy = proxy, quanti = quanti, ncores = ncores) return(exp) } #'@rdname DynBiasCorrection @@ -121,12 +120,11 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', #'dim (expL) <- c(time =100,lat = 4, lon = 5) #'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) #'dim (obsL) <- c(time = 100,lat = 4, lon = 5) -#'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, -#' proxy= "dim", quanti = 0.6) +#'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT', +#' proxy= "dim", quanti = 0.6,time_dim='time') #'@export DynBiasCorrection<- function(exp, obs, method = 'QUANT', - proxy = "dim", quanti, - time_dim = 'time', ncores = NULL){ + proxy = "dim", quanti, ncores = NULL){ if (is.null(obs)) { stop("Parameter 'obs' cannot be NULL.") } @@ -151,21 +149,62 @@ DynBiasCorrection<- function(exp, obs, method = 'QUANT', predyn.exp <- Predictability(dim = attractor.exp$dim, theta = attractor.exp$theta) + if (!(any(names(dim(exp)) %in% 'time'))){ + if (any(names(dim(exp)) %in% 'sdate')) { + if (any(names(dim(exp)) %in% 'ftime')) { + exp <- MergeDims(exp, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + } + if (!(any(names(dim(obs)) %in% 'time'))){ + if (any(names(dim(obs)) %in% 'sdate')) { + if (any(names(dim(obs)) %in% 'ftime')) { + obs <- MergeDims(obs, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + } + + + if(dim(obs)['member']!=dim(exp)['member']){ + names(dim(obs))[names(dim(obs))=='member'] <- 'memberObs' + } + + if(dim(obs)['dataset']!=dim(exp)['dataset']){ + names(dim(obs))[names(dim(obs))=='dataset'] <- 'datasetObs' + } if (proxy == "dim") { - adjusted <- Apply(list(exp, obs), target_dims = time_dim, + adjusted <- Apply(list(exp, obs), target_dims = 'time', fun = .dynbias, method, predyn.exp = predyn.exp$pred.dim$pos.d, predyn.obs = predyn.obs$pred.dim$pos.d, - ncores = ncores, output_dims = time_dim)$output1 + ncores = ncores, output_dims = 'time')$output1 } else if (proxy == "theta") { - adjusted <- Apply(list(exp, obs), target_dims = time_dim, + adjusted <- Apply(list(exp, obs), target_dims = 'time', fun = .dynbias, method, predyn.exp = predyn.exp$pred.theta$pos.t, predyn.obs = predyn.obs$pred.theta$pos.t, - ncores = ncores, output_dims = time_dim)$output1 + ncores = ncores, output_dims = 'time')$output1 } else { stop ("Parameter 'proxy' must be set as 'dim' or 'theta'.") } + + if(any(names(dim(adjusted)) %in% 'memberObs')){ + if(dim(adjusted)['memberObs'] == 1){ + adjusted <- Subset(adjusted,along='memberObs',indices=1,drop = 'selected') + }else{ + print('Dimension member in obs changed to memberObs') + } + } + + if(any(names(dim(adjusted)) %in% 'datasetObs')){ + if(dim(adjusted)['datasetObs'] == 1){ + adjusted <- Subset(adjusted,along='datasetObs',indices=1,drop = 'selected') + }else{ + print('Dimension dataset in obs changed to datasetObs') + } + } return(adjusted) } -- GitLab From 3eefddcb6bb313d7d082ff8d14f1c46e4090f507 Mon Sep 17 00:00:00 2001 From: nperez Date: Thu, 17 Jun 2021 17:50:47 +0200 Subject: [PATCH 15/19] updated docu and fixed common dimensions code --- R/CST_DynBiasCorrection.R | 34 ++++++++++++++++++++++------------ R/Predictability.R | 3 ++- man/CST_DynBiasCorrection.Rd | 5 +---- man/CST_ProxiesAttractor.Rd | 3 +-- man/DynBiasCorrection.Rd | 5 +---- man/Predictability.Rd | 6 ++++-- man/ProxiesAttractor.Rd | 20 +++++++++++--------- 7 files changed, 42 insertions(+), 34 deletions(-) diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R index 00d8dc68..52096ea6 100644 --- a/R/CST_DynBiasCorrection.R +++ b/R/CST_DynBiasCorrection.R @@ -31,7 +31,6 @@ #'"PTF","RQUANT","QUANT","SSPLIN" #'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. #'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta -#'@param time_dim a character string indicating the name of the temporal dimension #'@param ncores The number of cores to use in parallel computation #' #'@return dynbias an s2dvcube object with a bias correction performed @@ -56,7 +55,7 @@ #'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, #' Dates = list(start = time_obsL, end = time_obsL)) #'dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", -#' quanti = 0.6, time_dim = 'time') +#' quanti = 0.6) #' #'@export CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', @@ -107,7 +106,6 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', #'"PTF","RQUANT","QUANT","SSPLIN" #'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. #'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta -#'@param time_dim a character string indicating the name of the temporal dimension #'@param ncores The number of cores to use in parallel computation #' #'@return a multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta' @@ -121,7 +119,7 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', #'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) #'dim (obsL) <- c(time = 100,lat = 4, lon = 5) #'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT', -#' proxy= "dim", quanti = 0.6,time_dim='time') +#' proxy= "dim", quanti = 0.6) #'@export DynBiasCorrection<- function(exp, obs, method = 'QUANT', proxy = "dim", quanti, ncores = NULL){ @@ -165,15 +163,27 @@ DynBiasCorrection<- function(exp, obs, method = 'QUANT', } } } + + dim_exp <- dim(exp) + names_to_check <- names(dim_exp)[which(names(dim_exp) %in% + c('time', 'lat', 'lon', 'sdate') == FALSE)] + if (length(names_to_check) > 0) { + dim_obs <- dim(obs) + if (any(names(dim_obs) %in% names_to_check)) { + if (any(dim_obs[which(names(dim_obs) %in% names_to_check)] != + dim_exp[which(names(dim_exp) %in% names_to_check)])) { + for (i in names_to_check) { + pos <- which(names(dim_obs) == i) + names(dim(obs))[pos] <- ifelse(dim_obs[pos] != + dim_exp[which(names(dim_exp) == i)], + paste0('obs_', names(dim_obs[pos])), + names(dim(obs)[pos])) + } + warning("Common dimension names with different length are renamed.") + } + } + } - - if(dim(obs)['member']!=dim(exp)['member']){ - names(dim(obs))[names(dim(obs))=='member'] <- 'memberObs' - } - - if(dim(obs)['dataset']!=dim(exp)['dataset']){ - names(dim(obs))[names(dim(obs))=='dataset'] <- 'datasetObs' - } if (proxy == "dim") { adjusted <- Apply(list(exp, obs), target_dims = 'time', fun = .dynbias, method, diff --git a/R/Predictability.R b/R/Predictability.R index d17ea24c..73df65ee 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -25,6 +25,7 @@ #'for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints # #'@param theta list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes +#'@param ncores The number of cores to use in parallel computation #' #'@return A list of length 2: #' \itemize{ @@ -53,7 +54,7 @@ #'# imposing a threshold #' quanti <- 0.90 #'# computing dyn_scores from parameters dim and theta of the attractor -#' attractor <- ProxiesAttractor(dat = m, quanti = 0.60, iplot = FALSE) +#' attractor <- ProxiesAttractor(dat = m, quanti = 0.60) #' predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) #'@export #' diff --git a/man/CST_DynBiasCorrection.Rd b/man/CST_DynBiasCorrection.Rd index facf6f51..23eefd92 100644 --- a/man/CST_DynBiasCorrection.Rd +++ b/man/CST_DynBiasCorrection.Rd @@ -11,7 +11,6 @@ CST_DynBiasCorrection( method = "QUANT", proxy = "dim", quanti, - time_dim = "ftime", ncores = NULL ) } @@ -27,8 +26,6 @@ CST_DynBiasCorrection( \item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} -\item{time_dim}{a character string indicating the name of the temporal dimension} - \item{ncores}{The number of cores to use in parallel computation} } \value{ @@ -64,7 +61,7 @@ expL <- s2dv_cube(data = expL, lat = lat, lon = lon, obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, Dates = list(start = time_obsL, end = time_obsL)) dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", - quanti = 0.6, time_dim = 'time') + quanti = 0.6) } \references{ diff --git a/man/CST_ProxiesAttractor.Rd b/man/CST_ProxiesAttractor.Rd index 9757ff10..ecf9b9da 100644 --- a/man/CST_ProxiesAttractor.Rd +++ b/man/CST_ProxiesAttractor.Rd @@ -10,8 +10,7 @@ CST_ProxiesAttractor(data, quanti, ncores = NULL) \item{data}{a s2dv_cube object with the data to create the attractor. Must be a matrix with the timesteps in nrow and the grids in ncol(dat(time,grids)} -\item{quanti}{list of arbitrary length of secondary grids. Each secondary grid is to -be provided as a list of length 2 with longitudes and latitudes} +\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} \item{ncores}{The number of cores to use in parallel computation} } diff --git a/man/DynBiasCorrection.Rd b/man/DynBiasCorrection.Rd index bd60e0f3..f496c13d 100644 --- a/man/DynBiasCorrection.Rd +++ b/man/DynBiasCorrection.Rd @@ -11,7 +11,6 @@ DynBiasCorrection( method = "QUANT", proxy = "dim", quanti, - time_dim = "time", ncores = NULL ) } @@ -27,8 +26,6 @@ DynBiasCorrection( \item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} -\item{time_dim}{a character string indicating the name of the temporal dimension} - \item{ncores}{The number of cores to use in parallel computation} } \value{ @@ -50,7 +47,7 @@ expL <- rnorm(1:2000) dim (expL) <- c(time =100,lat = 4, lon = 5) obsL <- c(rnorm(1:1980),expL[1,,]*1.2) dim (obsL) <- c(time = 100,lat = 4, lon = 5) -dynbias <- DynBiasCorrection(exp = expL, obs = obsL, +dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT', proxy= "dim", quanti = 0.6) } \references{ diff --git a/man/Predictability.Rd b/man/Predictability.Rd index be0fb623..c36d6fa5 100644 --- a/man/Predictability.Rd +++ b/man/Predictability.Rd @@ -5,13 +5,15 @@ \title{Computing scores of predictability using two dinamical proxies based on dynamical systems theory.} \usage{ -Predictability(dim, theta) +Predictability(dim, theta, ncores = NULL) } \arguments{ \item{dim}{data to create the attractor. Must be a matrix with the timesteps in nrow and the grids in ncol for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints} \item{theta}{list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes} + +\item{ncores}{The number of cores to use in parallel computation} } \value{ A list of length 2: @@ -49,7 +51,7 @@ names(dim(m)) <- c('time', 'grid') # imposing a threshold quanti <- 0.90 # computing dyn_scores from parameters dim and theta of the attractor -attractor <- ProxiesAttractor(dat = m, quanti = 0.60, iplot = FALSE) +attractor <- ProxiesAttractor(dat = m, quanti = 0.60) predyn <- Predictability(dim = attractor$dim, theta = attractor$theta) } \references{ diff --git a/man/ProxiesAttractor.Rd b/man/ProxiesAttractor.Rd index b1fa71e3..768ba736 100644 --- a/man/ProxiesAttractor.Rd +++ b/man/ProxiesAttractor.Rd @@ -4,19 +4,12 @@ \alias{ProxiesAttractor} \title{Computing two dinamical proxies of the attractor.} \usage{ -ProxiesAttractor(data, quanti, iplot = FALSE, ncores = NULL) +ProxiesAttractor(data, quanti, ncores = NULL) } \arguments{ \item{data}{a multidimensional array with named dimensions to create the attractor. It requires a temporal dimension named 'time' and spatial dimensions called 'lat' and 'lon', or 'latitude' and 'longitude' or 'grid'.} -\item{quanti}{list of arbitrary length of secondary grids. Each secondary grid is to -be provided as a list of length 2 with longitudes and latitudes.} - -\item{iplot}{FALSE by default. If TRUE the function returns a simple plot -with three pannels on the top left a plot for local dimension 'd', on the top -right a plot for the inverse of the persistence 'theta', on the bottom a plot -for the 'attractor' plotted with the two properties 'd' (x axis) and 'theta' -(y axis)} +\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} \item{ncores}{The number of cores to use in parallel computation.} } @@ -37,6 +30,15 @@ Funtion based on the matlab code (davide.faranda@lsce.ipsl.fr) used in: mat <- array(rnorm(36 * 40), c(time = 36, grid = 40)) qm <- 0.90 # imposing a threshold Attractor <- ProxiesAttractor(data = mat, quanti = qm) +# to plot the result +time = c(1:length(Attractor$theta)) +layout(matrix(c(1, 3, 2, 3), 2, 2)) +plot(time, Attractor$dim, xlab = 'time', ylab = 'd', + main = 'local dimension', type = 'l') +plot(time, Attractor$theta, xlab = 'time', ylab = 'theta', main = 'theta') +plot(Attractor$dim, Attractor$theta, col = 'blue', + main = "Proxies of the Attractor", + xlab = "local dimension", ylab = "theta", lwd = 8, 'p') } \references{ -- GitLab From f3c9558813f253ca2b2a7a64e38e4335ff091af1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ver=C3=B3nica=20Torralba-Fern=C3=A1ndez?= Date: Fri, 25 Jun 2021 14:23:11 +0200 Subject: [PATCH 16/19] including a new parameter in CST_Dyn.. and checks in Predictability --- R/CST_DynBiasCorrection.R | 25 ++++++++++++++----------- R/Predictability.R | 21 ++++++++++++++++++++- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R index 52096ea6..bb7bc03b 100644 --- a/R/CST_DynBiasCorrection.R +++ b/R/CST_DynBiasCorrection.R @@ -29,6 +29,7 @@ #'@param obs an s2dv_cube object with the reference data #'@param method a character string indicating the method to apply bias correction among these ones: #'"PTF","RQUANT","QUANT","SSPLIN" +#'@param wetday logical indicating whether to perform wet day correction or not OR a numeric threshold below which all values are set to zero. #'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. #'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta #'@param ncores The number of cores to use in parallel computation @@ -58,7 +59,7 @@ #' quanti = 0.6) #' #'@export -CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', +CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=wetday, proxy = "dim", quanti, ncores = NULL) { if (!inherits(obs, 's2dv_cube')) { @@ -70,6 +71,7 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', "as output by CSTools::CST_Load.") } exp$data <- DynBiasCorrection(exp = exp$data, obs = obs$data, method = method, + wetday=wetday, proxy = proxy, quanti = quanti, ncores = ncores) return(exp) } @@ -104,6 +106,7 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', #'@param obs a multidimensional array with named dimensions with the observation data #'@param method a character string indicating the method to apply bias correction among these ones: #'"PTF","RQUANT","QUANT","SSPLIN" +#'@param wetday logical indicating whether to perform wet day correction or not OR a numeric threshold below which all values are set to zero. #'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. #'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta #'@param ncores The number of cores to use in parallel computation @@ -121,7 +124,7 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', #'dynbias <- DynBiasCorrection(exp = expL, obs = obsL, method='QUANT', #' proxy= "dim", quanti = 0.6) #'@export -DynBiasCorrection<- function(exp, obs, method = 'QUANT', +DynBiasCorrection<- function(exp, obs, method = 'QUANT',wetday=FALSE, proxy = "dim", quanti, ncores = NULL){ if (is.null(obs)) { stop("Parameter 'obs' cannot be NULL.") @@ -186,13 +189,13 @@ DynBiasCorrection<- function(exp, obs, method = 'QUANT', if (proxy == "dim") { adjusted <- Apply(list(exp, obs), target_dims = 'time', - fun = .dynbias, method, + fun = .dynbias, method, wetday, predyn.exp = predyn.exp$pred.dim$pos.d, predyn.obs = predyn.obs$pred.dim$pos.d, ncores = ncores, output_dims = 'time')$output1 } else if (proxy == "theta") { adjusted <- Apply(list(exp, obs), target_dims = 'time', - fun = .dynbias, method, + fun = .dynbias, method, wetday, predyn.exp = predyn.exp$pred.theta$pos.t, predyn.obs = predyn.obs$pred.theta$pos.t, ncores = ncores, output_dims = 'time')$output1 @@ -218,31 +221,31 @@ DynBiasCorrection<- function(exp, obs, method = 'QUANT', return(adjusted) } -.dynbias <- function(exp, obs, method, predyn.exp, predyn.obs) { +.dynbias <- function(exp, obs, method, wetday,predyn.exp, predyn.obs) { result <- array(rep(NA, length(exp))) res <- lapply(1:3, function(x) { exp_sub <- exp[predyn.exp[[x]]] obs_sub <- obs[predyn.obs[[x]]] - adjust <- .qbiascorrection(exp_sub, obs_sub, method) + adjust <- .qbiascorrection(exp_sub, obs_sub, method,wetday) result[predyn.exp[[x]]] <<- adjust return(NULL) }) return(result) } -.qbiascorrection <- function(expX, obsX, method) { +.qbiascorrection <- function(expX, obsX, method,wetday) { ## functions fitQmap and doQmap if (method == "PTF") { qm.fit <- fitQmap(obsX, expX, method = "PTF", transfun = "expasympt", - cost = "RSS", wett.day = TRUE) + cost = "RSS", wet.day = wetday) qmap <- doQmap(expX, qm.fit) } else if (method == "QUANT") { - qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01) + qm.fit <- fitQmap(obsX, expX, method = "QUANT", qstep = 0.01, wet.day = wetday) qmap <- doQmap(expX, qm.fit, type = "tricub") } else if (method == "RQUANT") { - qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01) + qm.fit <- fitQmap(obsX, expX, method = "RQUANT", qstep = 0.01,wet.day = wetday) qmap <- doQmap(expX, qm.fit, type = "linear") } else if (method == "SSPLIN") { - qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN") + qm.fit <- fitQmap(obsX, expX, qstep = 0.01, method = "SSPLIN",wet.day = wetday) qmap <- doQmap(expX, qm.fit) } else { stop ("Parameter 'method' doesn't match any of the available methods.") diff --git a/R/Predictability.R b/R/Predictability.R index 73df65ee..44d300a6 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -67,7 +67,26 @@ Predictability<- function(dim, theta, ncores = NULL) { # stop("Parameter 'theta' must be of the class 's2dv_cube', ", # "as output by CSTools::CST_Load.") # } - + if (any(names(dim(dim)) %in% 'sdate')) { + if (any(names(dim(dim)) %in% 'ftime')) { + dim <- MergeDims(dim, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(dim)) %in% 'time'))){ + stop("Parameter 'dim' must have a temporal dimension named 'time'.") + } + + if (any(names(dim(theta)) %in% 'sdate')) { + if (any(names(dim(theta)) %in% 'ftime')) { + theta <- MergeDims(theta, merge_dims = c('ftime', 'sdate'), + rename_dim = 'time') + } + } + if (!(any(names(dim(theta)) %in% 'time'))){ + stop("Parameter 'data' must have a temporal dimension named 'time'.") + } + pred <- Apply(list(dim, theta), target_dims = 'time', fun = .predictability, ncores = ncores) -- GitLab From 513dbcec21000b0f89e0ae77f7dc8d9f98e20704 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ver=C3=B3nica=20Torralba-Fern=C3=A1ndez?= Date: Fri, 25 Jun 2021 14:41:37 +0200 Subject: [PATCH 17/19] fixing typos --- R/Predictability.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/R/Predictability.R b/R/Predictability.R index 44d300a6..ddc516b8 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -1,5 +1,5 @@ #'@rdname Predictability -#'@title Computing scores of predictability using two dinamical proxies +#'@title Computing scores of predictability using two dynamical proxies #'based on dynamical systems theory. #' #'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it} @@ -11,7 +11,7 @@ #'computed with CST_ProxiesAttractor or ProxiesAttractor. These terciles will #'be used to measure the predictability of the system in dyn_scores. When the #'local dimension 'dim' is small and the inverse of persistence 'theta' is -#'small the predictability is higher, and viceversa. +#'small the predictability is high, and viceversa. #' #'@references Faranda, D., Alvarez-Castro, M.C., Messori, G., Rodriguez, D., #'and Yiou, P. (2019). The hammam effect or how a warm ocean enhances large @@ -21,7 +21,7 @@ #' Dynamical proxies of North Atlantic predictability and extremes. #' Scientific Reports, 7-41278, 2017. #' -#'@param dim data to create the attractor. Must be a matrix with the timesteps in nrow and the grids in ncol +#'@param dim data to create the attractor. must be a matrix with the timesteps in nrow and the grids in ncol #'for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints # #'@param theta list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes @@ -31,20 +31,20 @@ #' \itemize{ #' \item\code{pred.dim} {a list of two lists 'qdim' and 'pos.d'. The 'qdim' list #'contains values of local dimension 'dim' divided by terciles: -#'d1: lower tercile and more predictability, +#'d1: lower tercile (high predictability), #'d2: middle tercile, -#'d3: higher tercile and less predictability +#'d3: higher tercile (low predictability) #'The 'pos.d' list contains the position of each tercile in parameter 'dim'} #' #' \item\code{pred.theta} {a list of two lists 'qtheta' and 'pos.t'. #'The 'qtheta' list contains values of the inverse of persistence 'theta' #'divided by terciles: -#'th1: lower tercile and more predictability, +#'th1: lower tercile (high predictability), #'th2: middle tercile, -#'th3: higher tercile and less predictability +#'th3: higher tercile (low predictability) #'The 'pos.t' list contains the position of each tercile in parameter 'theta'} #'} -#'@return dyn_scores values from 0 to 1. A dyn_score of 1 indicates higher +#'@return dyn_scores values from 0 to 1. A dyn_score of 1 indicates the highest #'predictability. #' #'@examples -- GitLab From c4a575dfdd29e46f64b367e318b738ea5d422394 Mon Sep 17 00:00:00 2001 From: Carmen Alvarez-Castro Date: Fri, 2 Jul 2021 16:59:56 +0200 Subject: [PATCH 18/19] updating documentation --- R/CST_DynBiasCorrection.R | 73 ++++++++++++++++++++++++--------------- R/Predictability.R | 14 ++++---- 2 files changed, 53 insertions(+), 34 deletions(-) diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R index bb7bc03b..9aa33ced 100644 --- a/R/CST_DynBiasCorrection.R +++ b/R/CST_DynBiasCorrection.R @@ -27,11 +27,16 @@ #' #'@param exp an s2v_cube object with the experiment data #'@param obs an s2dv_cube object with the reference data -#'@param method a character string indicating the method to apply bias correction among these ones: -#'"PTF","RQUANT","QUANT","SSPLIN" -#'@param wetday logical indicating whether to perform wet day correction or not OR a numeric threshold below which all values are set to zero. -#'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. -#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta +#'@param method a character string indicating the method to apply bias +#'correction among these ones: "PTF","RQUANT","QUANT","SSPLIN" +#'@param wetday logical indicating whether to perform wet day correction +#'or not OR a numeric threshold below which all values are set to zero (by +#'default is set to 'FALSE'). +#'@param proxy a character string indicating the proxy for local dimension +#' 'dim' or inverse of persistence 'theta' to apply the dynamical +#' conditioned bias correction method. +#'@param quanti a number lower than 1 indicating the quantile to perform +#'the computation of local dimension and theta #'@param ncores The number of cores to use in parallel computation #' #'@return dynbias an s2dvcube object with a bias correction performed @@ -39,27 +44,30 @@ #' #'@examples #'# example 1: simple data s2dvcube style -#'set.seed(1) -#'expL <- rnorm(1:2000) -#'dim (expL) <- c(time =100,lat = 4, lon = 5) -#'obsL <- c(rnorm(1:1980),expL[1,,]*1.2) -#'dim (obsL) <- c(time = 100,lat = 4, lon = 5) -#'time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") -#'time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") -#'lon <- seq(-1,5,1.5) -#'lat <- seq(30,35,1.5) -# qm=0.98 # too high for this short dataset, it is possible that doesn't -#'# get the requirement, in that case it would be necessary select a lower qm -#'# for instance qm=0.60 -#'expL <- s2dv_cube(data = expL, lat = lat, lon = lon, +#' set.seed(1) +#' expL <- rnorm(1:2000) +#' dim (expL) <- c(time =100,lat = 4, lon = 5) +#' obsL <- c(rnorm(1:1980),expL[1,,]*1.2) +#' dim (obsL) <- c(time = 100,lat = 4, lon = 5) +#' time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") +#' time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") +#' lon <- seq(-1,5,1.5) +#' lat <- seq(30,35,1.5) +#' # qm=0.98 # too high for this short dataset, it is possible that doesn't# get the requirement, in that case it would be necessary select a lower qm +#' # for instance qm=0.60 +#' expL <- s2dv_cube(data = expL, lat = lat, lon = lon, #' Dates = list(start = time_expL, end = time_expL)) -#'obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, +#' obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, #' Dates = list(start = time_obsL, end = time_obsL)) -#'dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", +#' # to use DynBiasCorrection +#' dynbias1 <- DynBiasCorrection(exp = expL$data, obs = obsL$data, proxy= "dim", +#' quanti = 0.6) +#' # to use CST_DynBiasCorrection +#' dynbias2 <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", #' quanti = 0.6) #' #'@export -CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=wetday, +CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=FALSE, proxy = "dim", quanti, ncores = NULL) { if (!inherits(obs, 's2dv_cube')) { @@ -102,16 +110,25 @@ CST_DynBiasCorrection<- function(exp, obs, method = 'QUANT', wetday=wetday, #' Dynamical proxies of North Atlantic predictability and extremes. #' Scientific Reports, 7-41278, 2017. #' -#'@param exp a multidimensional array with named dimensions with the experiment data -#'@param obs a multidimensional array with named dimensions with the observation data -#'@param method a character string indicating the method to apply bias correction among these ones: +#'@param exp a multidimensional array with named dimensions with the +#'experiment data +#'@param obs a multidimensional array with named dimensions with the +#'observation data +#'@param method a character string indicating the method to apply bias +#'correction among these ones: #'"PTF","RQUANT","QUANT","SSPLIN" -#'@param wetday logical indicating whether to perform wet day correction or not OR a numeric threshold below which all values are set to zero. -#'@param proxy a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method. -#'@param quanti a number lower than 1 indicating the quantile to perform the computation of local dimension and theta +#'@param wetday logical indicating whether to perform wet day correction +#'or not OR a numeric threshold below which all values are set to zero (by +#'default is set to 'FALSE'). +#'@param proxy a character string indicating the proxy for local dimension +#''dim' or inverse of persistence 'theta' to apply the dynamical conditioned +#'bias correction method. +#'@param quanti a number lower than 1 indicating the quantile to perform the +#'computation of local dimension and theta #'@param ncores The number of cores to use in parallel computation #' -#'@return a multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta' +#'@return a multidimensional array with named dimensions with a bias correction +#'performed conditioned by local dimension 'dim' or inverse of persistence 'theta' #' #'@import multiApply #'@importFrom s2dverification Subset diff --git a/R/Predictability.R b/R/Predictability.R index ddc516b8..12cd6b41 100644 --- a/R/Predictability.R +++ b/R/Predictability.R @@ -21,10 +21,12 @@ #' Dynamical proxies of North Atlantic predictability and extremes. #' Scientific Reports, 7-41278, 2017. #' -#'@param dim data to create the attractor. must be a matrix with the timesteps in nrow and the grids in ncol -#'for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints +#'@param dim An array of N named dimensions containing the local dimension as +#'the output of CST_ProxiesAttractor or ProxiesAttractor. # -#'@param theta list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes +#'@param theta An array of N named dimensions containing the inverse of the +#'persistence 'theta' as the output of CST_ProxiesAttractor or ProxiesAttractor. +#' #'@param ncores The number of cores to use in parallel computation #' #'@return A list of length 2: @@ -92,9 +94,9 @@ Predictability<- function(dim, theta, ncores = NULL) { ncores = ncores) dim(pred$dyn_scores) <- dim(theta) return(list(pred.dim = list(qdim = list(pred$qdim.d1,pred$qdim.d2,pred$qdim.d3), - pos.d = list(pred$pos.d1,pred$pos.d2,pred$pos.d3)), - pred.theta = list(qtheta = list(pred$qtheta.th1,pred$qtheta.th2,pred$qtheta.th3), - pos.t = list(pred$pos.th1,pred$pos.th2,pred$pos.th3)), + pos.d = list(pred$pos.d1,pred$pos.d2,pred$pos.d3)), + pred.theta = list(qtheta = list(pred$qtheta.th1,pred$qtheta.th2,pred$qtheta.th3), + pos.t = list(pred$pos.th1,pred$pos.th2,pred$pos.th3)), dyn_scores = pred$dyn_scores)) } .predictability <- function(dim, theta) { -- GitLab From 548c902065b1626d364d3614211fa5f6ebf31a32 Mon Sep 17 00:00:00 2001 From: nperez Date: Mon, 13 Sep 2021 14:43:09 +0200 Subject: [PATCH 19/19] Tested on Win-builder --- DESCRIPTION | 2 +- NAMESPACE | 2 +- NEWS.md | 2 ++ R/CST_DynBiasCorrection.R | 3 ++- man/CST_DynBiasCorrection.Rd | 23 ++++++++++++++++++----- man/DynBiasCorrection.Rd | 24 ++++++++++++++++++------ man/Predictability.Rd | 21 +++++++++++---------- vignettes/Analogs_vignette.Rmd | 4 ++-- vignettes/RainFARM_vignette.Rmd | 2 +- 9 files changed, 56 insertions(+), 27 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index aecd89b3..97a80c72 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -84,4 +84,4 @@ VignetteBuilder: knitr License: Apache License 2.0 Encoding: UTF-8 LazyData: true -RoxygenNote: 7.0.2 +RoxygenNote: 7.0.1 diff --git a/NAMESPACE b/NAMESPACE index 26042662..34851f66 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -120,4 +120,4 @@ importFrom(utils,head) importFrom(utils,read.table) importFrom(utils,tail) importFrom(verification,verify) -useDynLib(CSTools) +useDynLib(CSTools, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index 9a4c78e0..154c6307 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,8 @@ **Submission date to CRAN: XX-06-2021** - New features: + + Dynamical Bias Correction method: `CST_ProxiesAttractors` and `CST_DynBiasCorrection` + (optionally `Predictability`) + CST_BiasCorrection and BiasCorrection allows to calibrate a forecast given the calibration in the hindcast by using parameter 'exp_cor'. - Fixes: diff --git a/R/CST_DynBiasCorrection.R b/R/CST_DynBiasCorrection.R index 9aa33ced..20c263c6 100644 --- a/R/CST_DynBiasCorrection.R +++ b/R/CST_DynBiasCorrection.R @@ -53,7 +53,8 @@ #' time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") #' lon <- seq(-1,5,1.5) #' lat <- seq(30,35,1.5) -#' # qm=0.98 # too high for this short dataset, it is possible that doesn't# get the requirement, in that case it would be necessary select a lower qm +#' # qm=0.98 # too high for this short dataset, it is possible that doesn't +#' # get the requirement, in that case it would be necessary select a lower qm #' # for instance qm=0.60 #' expL <- s2dv_cube(data = expL, lat = lat, lon = lon, #' Dates = list(start = time_expL, end = time_expL)) diff --git a/man/CST_DynBiasCorrection.Rd b/man/CST_DynBiasCorrection.Rd index 23eefd92..e2467e72 100644 --- a/man/CST_DynBiasCorrection.Rd +++ b/man/CST_DynBiasCorrection.Rd @@ -9,6 +9,7 @@ CST_DynBiasCorrection( exp, obs, method = "QUANT", + wetday = FALSE, proxy = "dim", quanti, ncores = NULL @@ -19,12 +20,19 @@ CST_DynBiasCorrection( \item{obs}{an s2dv_cube object with the reference data} -\item{method}{a character string indicating the method to apply bias correction among these ones: -"PTF","RQUANT","QUANT","SSPLIN"} +\item{method}{a character string indicating the method to apply bias +correction among these ones: "PTF","RQUANT","QUANT","SSPLIN"} -\item{proxy}{a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method.} +\item{wetday}{logical indicating whether to perform wet day correction +or not OR a numeric threshold below which all values are set to zero (by +default is set to 'FALSE').} -\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} +\item{proxy}{a character string indicating the proxy for local dimension +'dim' or inverse of persistence 'theta' to apply the dynamical +conditioned bias correction method.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform +the computation of local dimension and theta} \item{ncores}{The number of cores to use in parallel computation} } @@ -54,13 +62,18 @@ time_obsL <- paste(rep("01", 100), rep("01", 100), 1920 : 2019, sep = "-") time_expL <- paste(rep("01", 100), rep("01", 100), 1929 : 2019, sep = "-") lon <- seq(-1,5,1.5) lat <- seq(30,35,1.5) +# qm=0.98 # too high for this short dataset, it is possible that doesn't # get the requirement, in that case it would be necessary select a lower qm # for instance qm=0.60 expL <- s2dv_cube(data = expL, lat = lat, lon = lon, Dates = list(start = time_expL, end = time_expL)) obsL <- s2dv_cube(data = obsL, lat = lat, lon = lon, Dates = list(start = time_obsL, end = time_obsL)) -dynbias <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", +# to use DynBiasCorrection +dynbias1 <- DynBiasCorrection(exp = expL$data, obs = obsL$data, proxy= "dim", + quanti = 0.6) +# to use CST_DynBiasCorrection +dynbias2 <- CST_DynBiasCorrection(exp = expL, obs = obsL, proxy= "dim", quanti = 0.6) } diff --git a/man/DynBiasCorrection.Rd b/man/DynBiasCorrection.Rd index f496c13d..e6de373c 100644 --- a/man/DynBiasCorrection.Rd +++ b/man/DynBiasCorrection.Rd @@ -9,27 +9,39 @@ DynBiasCorrection( exp, obs, method = "QUANT", + wetday = FALSE, proxy = "dim", quanti, ncores = NULL ) } \arguments{ -\item{exp}{a multidimensional array with named dimensions with the experiment data} +\item{exp}{a multidimensional array with named dimensions with the +experiment data} -\item{obs}{a multidimensional array with named dimensions with the observation data} +\item{obs}{a multidimensional array with named dimensions with the +observation data} -\item{method}{a character string indicating the method to apply bias correction among these ones: +\item{method}{a character string indicating the method to apply bias +correction among these ones: "PTF","RQUANT","QUANT","SSPLIN"} -\item{proxy}{a character string indicating the proxy for local dimension 'dim' or inverse of persistence 'theta' to apply the dynamical conditioned bias correction method.} +\item{wetday}{logical indicating whether to perform wet day correction +or not OR a numeric threshold below which all values are set to zero (by +default is set to 'FALSE').} -\item{quanti}{a number lower than 1 indicating the quantile to perform the computation of local dimension and theta} +\item{proxy}{a character string indicating the proxy for local dimension +'dim' or inverse of persistence 'theta' to apply the dynamical conditioned +bias correction method.} + +\item{quanti}{a number lower than 1 indicating the quantile to perform the +computation of local dimension and theta} \item{ncores}{The number of cores to use in parallel computation} } \value{ -a multidimensional array with named dimensions with a bias correction performed conditioned by local dimension 'dim' or inverse of persistence 'theta' +a multidimensional array with named dimensions with a bias correction +performed conditioned by local dimension 'dim' or inverse of persistence 'theta' } \description{ This function perform a bias correction conditioned by the diff --git a/man/Predictability.Rd b/man/Predictability.Rd index c36d6fa5..d37efcdc 100644 --- a/man/Predictability.Rd +++ b/man/Predictability.Rd @@ -2,16 +2,17 @@ % Please edit documentation in R/Predictability.R \name{Predictability} \alias{Predictability} -\title{Computing scores of predictability using two dinamical proxies +\title{Computing scores of predictability using two dynamical proxies based on dynamical systems theory.} \usage{ Predictability(dim, theta, ncores = NULL) } \arguments{ -\item{dim}{data to create the attractor. Must be a matrix with the timesteps in nrow and the grids in ncol -for instance: dat(time,grids)=(1:24418,1:1060), where time=24418 timesteps and grids=1060 gridpoints} +\item{dim}{An array of N named dimensions containing the local dimension as +the output of CST_ProxiesAttractor or ProxiesAttractor.} -\item{theta}{list of arbitrary length of secondary grids. Each secondary grid is to be provided as a list of length 2 with longitudes and latitudes} +\item{theta}{An array of N named dimensions containing the inverse of the +persistence 'theta' as the output of CST_ProxiesAttractor or ProxiesAttractor.} \item{ncores}{The number of cores to use in parallel computation} } @@ -20,21 +21,21 @@ A list of length 2: \itemize{ \item\code{pred.dim} {a list of two lists 'qdim' and 'pos.d'. The 'qdim' list contains values of local dimension 'dim' divided by terciles: -d1: lower tercile and more predictability, +d1: lower tercile (high predictability), d2: middle tercile, -d3: higher tercile and less predictability +d3: higher tercile (low predictability) The 'pos.d' list contains the position of each tercile in parameter 'dim'} \item\code{pred.theta} {a list of two lists 'qtheta' and 'pos.t'. The 'qtheta' list contains values of the inverse of persistence 'theta' divided by terciles: -th1: lower tercile and more predictability, +th1: lower tercile (high predictability), th2: middle tercile, -th3: higher tercile and less predictability +th3: higher tercile (low predictability) The 'pos.t' list contains the position of each tercile in parameter 'theta'} } -dyn_scores values from 0 to 1. A dyn_score of 1 indicates higher +dyn_scores values from 0 to 1. A dyn_score of 1 indicates the highest predictability. } \description{ @@ -42,7 +43,7 @@ This function divides in terciles the two dynamical proxies computed with CST_ProxiesAttractor or ProxiesAttractor. These terciles will be used to measure the predictability of the system in dyn_scores. When the local dimension 'dim' is small and the inverse of persistence 'theta' is -small the predictability is higher, and viceversa. +small the predictability is high, and viceversa. } \examples{ # Creating an example of matrix dat(time,grids): diff --git a/vignettes/Analogs_vignette.Rmd b/vignettes/Analogs_vignette.Rmd index 5a52a05d..31fb9c29 100644 --- a/vignettes/Analogs_vignette.Rmd +++ b/vignettes/Analogs_vignette.Rmd @@ -156,7 +156,7 @@ dateseq <- format(seq(start, end, by = "year"), "%Y%m%d") Using the `CST_Load` function from **CSTool package**, the data available in our data store can be loaded. The following lines show how this function can be used. The experimental datasets are interpolated to the ERA5 grid by specifying the 'grid' parameter while ERA5 doesn't need to be interpolated. While parameter leadtimemax is set to 1 for the experimental dataset, it is set to 31 for the observations, returning the daily observations for October for the years requested in 'sdate' (2000-2006). -Download the data to run the recipe in the link https://downloads.cmcc.bo.it/d_chaves/ANALOGS/data_for_Analogs.Rdat or ask carmen.alvarez-castro at cmcc.it or nuria.perez at bsc.es. +Download the data to run the recipe under the HTTPS: downloads.cmcc.bo.it/d_chaves/ANALOGS/data_for_Analogs.Rdat or ask carmen.alvarez-castro at cmcc.it or nuria.perez at bsc.es. ``` exp <- list(name = 'ECMWF_system4_m1', @@ -373,4 +373,4 @@ down4 <- CST_Analogs(expL = expPSL, obsL = obsPSL, AnalogsInfo = TRUE, In this case, the best analog is still being 7th of October, 2005. -*Note: You can compute the anomalies values before applying the criterias (as in Yiou et al, 2013) using `CST_Anomaly` of CSTools package* \ No newline at end of file +*Note: You can compute the anomalies values before applying the criterias (as in Yiou et al, 2013) using `CST_Anomaly` of CSTools package* diff --git a/vignettes/RainFARM_vignette.Rmd b/vignettes/RainFARM_vignette.Rmd index dbcb48a4..c47d0e73 100644 --- a/vignettes/RainFARM_vignette.Rmd +++ b/vignettes/RainFARM_vignette.Rmd @@ -118,7 +118,7 @@ RainFARM has downscaled the original field with a realistic fine-scale correlati The area of interest in our example presents a complex orography, but the basic RainFARM algorithm used does not consider topographic elevation in deciding how to distribute fine-scale precipitation. A long term climatology of the downscaled fields would have a resolution comparable to that of the original coarse fields and would not resemble the fine-scale structure of an observed climatology. If an external fine-scale climatology of precipitation is available, we can use the method discussed in Terzago et al. (2018) to change the distribution of precipitation by RainFARM for each timestep, so that the long-term average is close to this reference climatology in terms of precipitation distribution (while the total precipitation amount of the original fields to downscale is preserved). -Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (https://www.worldclim.org) or CHELSA (https://chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://gdal.org). +Suitable climatology files could be for example a fine-scale precipitation climatology from a high-resolution regional climate model (see e.g. Terzago et al. 2018), a local high-resolution gridded climatology from observations, or a reconstruction such as those which can be downloaded from the WORLDCLIM (https://www.worldclim.org) or CHELSA (chelsa-climate.org) websites. The latter data will need to be converted to NetCDF format before being used (see for example the GDAL tools (https://gdal.org). We will assume that a copy of the WORLDCLIM precipitation climatology at 30 arcseconds (about 1km resolution) is available in the local file `medscope.nc`. From this file we can derive suitable weights to be used with RainFARM using the `CST_RFWeights` functions as follows: ```{r} ww <- CST_RFWeights("./worldclim.nc", nf = 20, lon = exp$lon, lat = exp$lat) -- GitLab