Newer
Older
EVA RIFA ROVIRA
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
# WIP
PeriodStandardization <- function(data, data_cor = NULL, time_dim = 'syear',
leadtime_dim = 'time', memb_dim = 'ensemble',
ref_period = NULL, cross_validation = FALSE,
handle_infinity = FALSE, param_error = -9999,
method = 'parametric', distribution = 'log-Logistic',
na.rm = FALSE, ncores = NULL) {
# Initial checks
target_dims <- c(leadtime_dim, time_dim, memb_dim)
if (is.null(ref_period)) {
ref.start <- NULL
ref.end <- NULL
} else {
ref.start <- ref_period[[1]]
ref.end <- ref_period[[2]]
}
# Standardization
if (is.null(data_cor)) {
spei <- Apply(data = list(data),
target_dims = target_dims,
fun = .standardization,
leadtime_dim = leadtime_dim,
time_dim = time_dim, memb_dim = memb_dim,
ref_period = ref_period, handle_infinity = handle_infinity,
cross_validation = cross_validation, param_error = param_error,
method = method, distribution = distribution,
na.rm = na.rm, ncores = ncores)$output1
} else {
spei <- Apply(data = list(data, data_cor), target_dims = target_dims,
fun = .standardization,
leadtime_dim = leadtime_dim,
time_dim = time_dim, memb_dim = memb_dim,
ref_period = ref_period, handle_infinity = handle_infinity,
cross_validation = cross_validation, param_error = param_error,
method = method, distribution = distribution,
na.rm = na.rm, ncores = ncores)$output1
}
return(spei)
}
# data <- array(rnorm(10), c(time = 3, syear = 6, ensemble = 25))
# res <- .standardization(data = data)
# data_subset <- ClimProjDiags::Subset(data, along = leadtime_dim,
# indices = ff, drop = 'selected')
.standardization <- function(data, data_cor = NULL, leadtime_dim = 'time',
time_dim = 'syear', memb_dim = 'ensemble',
ref_period = NULL, handle_infinity = FALSE,
cross_validation = FALSE, param_error = -9999,
method = 'parametric', distribution = 'log-Logistic',
na.rm = FALSE) {
# data: [leadtime_dim, time_dim, memb_dim]
# maximum number of parameters needed
nleadtime <- as.numeric(dim(data)[leadtime_dim])
ntime <- as.numeric(dim(data)[time_dim])
nmemb <- as.numeric(dim(data)[memb_dim])
fit = 'ub-pwm'
coef = switch(distribution,
"Gamma" = array(NA, dim = 2, dimnames = list(c('alpha', 'beta'))),
"log-Logistic" = array(NA, dim = 3, dimnames = list(c('xi', 'alpha', 'kappa'))),
"PearsonIII" = array(NA, dim = 3, dimnames = list(c('mu', 'sigma', 'gamma'))))
if (is.null(data_cor)) {
# cross_val = TRUE
spei_mod <- data*NA
print(dim(spei_mod))
for (ff in 1:dim(data)[leadtime_dim]) {
data2 <- data[ff, , ]
params_result <- array(dim = c(dim(data)[time_dim], length(coef)))
print(dim(data2))
if (!is.null(ref.start) && !is.null(ref.end)) {
data.fit <- window(data2, ref.start, ref.end)
} else {
data.fit <- data2
}
for (nsd in 1:dim(data)[time_dim]) {
acu <- as.vector(data.fit[-nsd, ])
acu.sorted <- sort.default(acu, method = "quick")
if (na.rm) {
acu.sorted <- acu.sorted[!is.na(acu.sorted)]
}
if (!any(is.na(acu.sorted))) {
print('Inside acu.sorted')
if (length(acu.sorted) != 0) {
acu_sd <- sd(acu.sorted)
if (!is.na(acu_sd) & acu_sd != 0) {
if (distribution != "log-Logistic") {
pze <- sum(acu == 0) / length(acu)
acu.sorted <- acu.sorted[acu.sorted > 0]
}
if (length(acu.sorted) >= 4) {
print('acu.sorted')
print(acu.sorted)
pwm = switch(fit,
"pp-pwm" = pwm.pp(acu.sorted, -0.35, 0, nmom = 3),
pwm.ub(acu.sorted, nmom = 3)
# TLMoments::PWM(acu.sorted, order = 0:2)
)
lmom <- pwm2lmom(pwm)
if (!(!are.lmom.valid(lmom) || anyNA(lmom[[1]]) || any(is.nan(lmom[[1]])))) {
fortran_vec = c(lmom$lambdas[1:2], lmom$ratios[3])
f_params = switch(distribution,
"log-Logistic" = tryCatch(lmom::pelglo(fortran_vec),
error = function(e){parglo(lmom)$para}),
"Gamma" = tryCatch(lmom::pelgam(fortran_vec),
error = function(e){pargam(lmom)$para}),
"PearsonIII" = tryCatch(lmom::pelpe3(fortran_vec),
error = function(e){parpe3(lmom)$para}))
if (distribution == 'log-Logistic' && fit == 'max-lik') {
f_params = parglo.maxlik(acu.sorted, f_params)$para
}
params_result[nsd, ] <- f_params
}
}
if (all(is.na(params_result[nsd,]))) {
cdf_res <- NA
} else {
f_params <- params_result[nsd,]
f_params <- f_params[which(!is.na(f_params))]
cdf_res = switch(distribution,
"log-Logistic" = lmom::cdfglo(data, f_params),
"Gamma" = lmom::cdfgam(data, f_params),
"PearsonIII" = lmom::cdfpe3(data, f_params))
}
std_index_cv <- array(qnorm(cdf_res), dim = c(ntime, nmemb))
spei_mod[ff, nsd, ] <- std_index_cv[nsd, ]
}
}
}
}
}
} else {
# cross_val = FALSE
}
return(spei_mod)
}