Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#'@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<qd3)
qdim <- list(d1=dim[d1],d2=dim[d2],d3=dim[d3])
pos.d <- list(d1=d1,d2=d2,d3=d3)
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#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)
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)]
d1th2 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th2)]
d2th1 <- pos.d$d2[which(pos.d$d2 %in% pos.t$th1)]
d1th3 <- pos.d$d1[which(pos.d$d1 %in% pos.t$th3)]
d3th1 <- pos.d$d3[which(pos.d$d3 %in% pos.t$th1)]
d2th2 <- pos.d$d2[which(pos.d$d2 %in% pos.t$th2)]
d2th3 <- pos.d$d2[which(pos.d$d2 %in% pos.t$th3)]
d3th2 <- pos.d$d3[which(pos.d$d3 %in% pos.t$th2)]
d3th3 <- pos.d$d3[which(pos.d$d3 %in% pos.t$th3)]
#scores values
dyn_scores <- c(1:length(dim))
dyn_scores[which(pos %in% d1th1)]<- 9/9
dyn_scores[which(pos %in% d1th2)]<- 8/9
dyn_scores[which(pos %in% d2th1)]<- 7/9
dyn_scores[which(pos %in% d1th3)]<- 6/9
dyn_scores[which(pos %in% d3th1)]<- 5/9
dyn_scores[which(pos %in% d2th2)]<- 4/9
dyn_scores[which(pos %in% d2th3)]<- 3/9
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))
}