# R/decomps.R In hanel/muscabicor:

#### Defines functions run_gr4jrun_bilgr4jplot_gbias2plot_gbiasplot_bias_smryplot_biasmelt_bias2melt_biasgbias2gbiasbiasget_valuegive_mapdecomprev_difdifaggprob

```prob = function(x){
(rank(x)-.3) / (length(x)+.4)
}

agg = function(x, variable, level = '?',  ...){
fun = if (variable=='PR' & level!='G1..') (sum) else (mean)
fun(x, ...)
}

dif = function(x, y, variable){
if (variable == 'TAS') return(x-y)
r = x/y
r[y==0] = 0
return(r)
}

rev_dif = function(x, y, variable){
if (variable == 'TAS') return(x+y)
r = x*y
#r[y==0] = 0
return(r)
}

decomp = function(dta, f = formula(~ Y + M + D), year_starts = -months(2), complete = TRUE, wdt = 0.1){
dta = copy(dta)
m = give_map(dta[, .(DTM)], f = f, year_starts = year_starts)
mdta = melt(dta[m, on = 'DTM'], measure = names(dta)[names(dta)!='DTM'])

r = mdta[, .(DTM, variable, value)]
r[variable=='PR' & value < wdt, value := 0]

for (i in which(names(mdta) %in% attr(m, 'terms')\$pterms)){
n = names(mdta)[i]
if (grepl('Y|G', n)) p = mdta[,  .(DTM, agg(value, variable[1], level = n)), by = .(variable, eval(parse(text = n)))]
if (grepl('M', n)) p = mdta[,  .(DTM, agg(value, variable[1], level = n)), by = .(variable, Y1, eval(parse(text = n)))]
if (grepl('D.', n)) p = mdta[,  .(DTM, agg(value, variable[1], level = n)), by = .(variable, Y1, M1, eval(parse(text = n)))]
if (!grepl('D.', n) & n!='D1' & grepl('D', n)) p = mdta[,  .(DTM, agg(value, variable[1], level = n)), by = .(variable, eval(parse(text = n)))]
if (n=='D1') p = mdta[, .(DTM, V2 = value)]
r[, (tolower(n)) := p\$V2]

}
r = m[r, on = 'DTM']

for (i in which(names(r) %in% tolower(attr(m, 'terms')\$pterms))[-1]){
r[, (paste0('w_', names(r)[i])) := .(dif(eval(parse(text = names(r)[i])), eval(parse(text = names(r)[i-1])), variable =  variable[1])), by = variable]
}

if (complete == TRUE){
from = r[1, DTM] %m+% (months(12) + year_starts)
to = r[nrow(r), DTM] %m+% year_starts
r = r[(DTM >= from) & (DTM <= (to))]
}
structure(copy(r), terms = attr(m, 'terms'), formula = f)

}

give_map = function(map, f , year_starts){

#map = copy(map)

terms = c('G*1', gsub(' ', '', strsplit(as.character(f)[2], '\\+')[[1]]))
scale = gsub('(\\d+)|\\*|\\/', '', terms)
len = suppressWarnings(as.integer(gsub("[^\\d]+", "", terms, perl = TRUE)))
len[is.na(len)] = 1
fce = gsub('[^[:punct:],]', '', terms)
fce[fce==''] = '*'

terms = paste0(len, fce, scale)
breaks = gsub('\\*', ' ', terms)
breaks = gsub('Y', 'year', breaks)
breaks = gsub('M', 'month', breaks)
breaks = gsub('D', 'day', breaks)
pterms = paste0(scale, gsub('\\*', '', fce), len)
pterms = gsub('/', '.', pterms)

i = 1
map[, G1:=1]

for (i in which(scale=='Y')){
map[, c(pterms[i]) :=  as.integer(cut(DTM %m-% year_starts, breaks = breaks[i]))]
}

for (i in which(scale=='M')){
map[, c(pterms[i]) :=  month(cut(DTM %m-% year_starts, breaks = breaks[i])), by = Y1]
}

for (i in which(scale=='D' & fce=='*' & len!=1)){
map[, c(pterms[i]) :=  1 + (days(DTM) - days(DTM)[1]) %/% days(len[i])]
}

for (i in which(scale=='D' & fce=='/' & len!=1)){
map[, c(pterms[i]) := as.integer(cut(DTM, breaks = len[i])), by = .(Y1, M1)]
}

map[, D1:= 1 + as.integer(DTM) - as.integer(DTM[1])]

structure(copy(map), terms = list(terms = terms, pterms = pterms, scale = scale, fce = fce, len = len))
}

get_value = function(d, level = 'Y1', daily_by = 'M3', w = FALSE){
l = if (w==FALSE) {level} else {paste0('w_', level)}
if (grepl('Y|G', level)) r = d[, eval(parse(text = tolower(l)))[1], by = .(variable, eval(parse(text = toupper(level)))) ]
if (grepl('M', level)) r = d[, eval(parse(text = tolower(l)))[1], by = .(variable, Y1, eval(parse(text = toupper(level)))) ]
if (grepl('D.', level)) r = d[, eval(parse(text = tolower(l)))[1], by = .(variable, Y1, db = eval(parse(text = daily_by)), eval(parse(text = toupper(level)))) ]
if (!grepl('D.', level) & level!='D1' & grepl('D', level)) r = d[, eval(parse(text = tolower(l)))[1], by = .(variable, eval(parse(text = toupper(level)))) ]
if (level=='D1') r = d[, eval(parse(text = tolower(l)))[1], by = .(variable, db = eval(parse(text = daily_by)), eval(parse(text = toupper(level)))) ]
if ('db' %in% names(r)) setnames(r, 'db', daily_by)
setnames(r, c('parse', 'V1'), c(level, tolower(l)))
copy(r)
}

bias = function(obs, ctrl, w = FALSE, daily_by = 'M1', wet_int_only = TRUE, wdt = 0.1){

pt = attr(obs, 'terms')\$pterms
if (w) pt = pt[pt!='G1']
#r = obs[, .(DTM, variable, value)]
i = pt[3]

for (i in pt){

p = seq(.05, .95, by = .05)#if (!grepl('M1|D', i)) (seq(.1, .9, by = .1)) else (seq(.1, .9, by = .05))
xo = get_value(obs, i, daily_by = daily_by, w = w)
xc = get_value(ctrl, i, daily_by = daily_by, w = w)
ii = if (!w) (tolower(i)) else (paste0('w_', tolower(i)))

if (grepl('Y|G', i)){

qxo = xo[, .(p, O = quantile(eval(parse(text = ii )), p = p)), by = .(variable)]
qxc = xc[, .(p, C = quantile(eval(parse(text = ii )), p = p)), by = .(variable)]
b = qxo[qxc, on = c('variable', 'p')]
}

if (grepl('M', i)){
qxo = xo[, .(p, O = quantile(eval(parse(text = ii )), p = p)), by = .(variable, eval(parse(text = i)))]
qxc = xc[, .(p, C = quantile(eval(parse(text = ii )), p = p)), by = .(variable, eval(parse(text = i)))]
b = qxo[qxc, on = c('variable', 'parse', 'p')]
}

if (grepl('D', i) & !grepl('D1', i)){
iii = if (wet_int_only) {paste0(ii, '[', ii, '>', wdt, ']')} else {ii}
qxo = xo[, .(p, O = quantile(eval(parse(text = iii )), p = p)), by = .(variable, db = eval(parse(text = daily_by)))]
qxc = xc[, .(p, C = quantile(eval(parse(text = iii )), p = p)), by = .(variable, db = eval(parse(text = daily_by)))]
b = qxo[qxc, on = c('variable', 'db',  'p')]
setnames(b, 'db', daily_by)
}

if (grepl('D1', i)){
iii = if (wet_int_only) {paste0(ii, '[', ii, '>', wdt, ']')} else {ii}
#if (!w){
qxo = xo[, .(p, O = quantile(eval(parse(text = iii)), p =p)), by = .(variable, db = eval(parse(text = daily_by)))]
qxc = xc[, .(p, C = quantile(eval(parse(text = iii)), p = p)), by = .(variable, db = eval(parse(text = daily_by)))]
# #} else {
#   qxo = xo[, .(p, O = quantile(w_d1, p =p)), by = .(variable, db = eval(parse(text = daily_by)))]
#   qxc = xc[, .(p, C = quantile(w_d1, p = p)), by = .(variable, db = eval(parse(text = daily_by)))]
# }
b = qxo[qxc, on = c('variable', 'db', 'p')]
setnames(b, 'db', daily_by)
}

b[variable != 'TAS', dif := C/O]
b[variable == 'TAS', dif := C-O]
if ('parse' %in% names(b)) setnames(b, 'parse', i)

if ( (i=='G1'&w==FALSE)|(w==TRUE & i==pt[1]) ) {
ppt = unique(c(pt[grepl('G|M', pt)], daily_by))
ppt = c('variable', ppt)#, pt[!(pt %in% ppt) & !(grepl('D1|Y', pt ))])
RES = obs[, .(p), by = ppt]
RES = b[, .(variable, p, dif)][RES, on = c('variable', 'p')]
setnames(RES, 'dif', if(!w) ('g1') else ('w_y1') )
} else {

b[, c('O', 'C') := NULL]
key = names(b)[names(b) %in% names(RES)]
RES = b[RES, on = key]
setnames(RES, 'dif', ii)

}

}

setcolorder(RES, c('variable', pt[!grepl('Y|D', pt)], c('p', paste0(if (w) {'w_'} else {''}, tolower(pt)))))
copy(RES)
}

tscale = Vectorize(function(x, nyears = 30){

if (x=='g1') return(as.double(duration(nyears, 'year'))/3600)
x = gsub('w_', '',  x)
x = gsub('y', 'year', x)
x = gsub('m', 'month', x)
x = gsub('d', 'day', x)
num = suppressWarnings(as.integer(gsub("[^\\d]+", "", x, perl = TRUE)))
if (grepl('\\.', x)) {num = 1/num; x = 'month'}
#x = gsub('\\.', '', x)
x = gsub('(\\d+)|\\.', '', x)
as.double(duration(num, x))/3600

})

gbias = function(obs, ctrl, w = FALSE, daily_by = 'M1', wet_int_only = TRUE, p = seq(.1, .9, by = .1),  fun = function(x)quantile(x, p = p), wdt = 0.1){

pt = attr(obs, 'terms')\$pterms
if (w) pt = pt[pt!='G1']
#r = obs[, .(DTM, variable, value)]
i = pt[3]

for (i in pt){

#p = seq(.05, .95, by = .05)#if (!grepl('M1|D', i)) (seq(.1, .9, by = .1)) else (seq(.1, .9, by = .05))
xo = get_value(obs, i, daily_by = daily_by, w = w)
xc = get_value(ctrl, i, daily_by = daily_by, w = w)
ii = if (!w) (tolower(i)) else (paste0('w_', tolower(i)))

if (grepl('Y|G', i)){

qxo = xo[, .(p, O = fun(eval(parse(text = ii )))), by = .(variable)]
qxc = xc[, .(p, C = fun(eval(parse(text = ii )))), by = .(variable)]
b = qxo[qxc, on = c('variable', 'p')]
}

if (grepl('M', i)){
qxo = xo[, .(p, O = fun(eval(parse(text = ii )))), by = .(variable, eval(parse(text = i)))]
qxc = xc[, .(p, C = fun(eval(parse(text = ii )))), by = .(variable, eval(parse(text = i)))]
b = qxo[qxc, on = c('variable', 'parse', 'p')]
}

if (grepl('D', i) & !grepl('D1', i)){
iii = if (wet_int_only) {paste0(ii, '[', ii, '>', wdt, ']')} else {ii}
qxo = xo[, .(p, O = fun(eval(parse(text = iii )))), by = .(variable, db = eval(parse(text = daily_by)))]
qxc = xc[, .(p, C = fun(eval(parse(text = iii )))), by = .(variable, db = eval(parse(text = daily_by)))]
b = qxo[qxc, on = c('variable', 'db',  'p')]
setnames(b, 'db', daily_by)
}

if (grepl('D1', i)){
iii = if (wet_int_only) {paste0(ii, '[', ii, '>', wdt, ']')} else {ii}
#if (!w){
qxo = xo[, .(p, O = fun(eval(parse(text = iii)))), by = .(variable, db = eval(parse(text = daily_by)))]
qxc = xc[, .(p, C = fun(eval(parse(text = iii)))), by = .(variable, db = eval(parse(text = daily_by)))]
# #} else {
#   qxo = xo[, .(p, O = quantile(w_d1, p =p)), by = .(variable, db = eval(parse(text = daily_by)))]
#   qxc = xc[, .(p, C = quantile(w_d1, p = p)), by = .(variable, db = eval(parse(text = daily_by)))]
# }
b = qxo[qxc, on = c('variable', 'db', 'p')]
setnames(b, 'db', daily_by)
}

b[! variable %in% c('TAS', 'Xd'), dif := C/O]
b[variable %in% c('TAS', 'Xd'), dif := C-O]
if ('parse' %in% names(b)) setnames(b, 'parse', i)

if ( (i=='G1'&w==FALSE)|(w==TRUE & i==pt[1]) ) {
ppt = unique(c(pt[grepl('G|M', pt)], daily_by))
ppt = c('variable', ppt)#, pt[!(pt %in% ppt) & !(grepl('D1|Y', pt ))])
RES = obs[, .(p), by = ppt]
RES = b[, .(variable, p, dif)][RES, on = c('variable', 'p')]
setnames(RES, 'dif', if(!w) ('g1') else ('w_y1') )
} else {

b[, c('O', 'C') := NULL]
key = names(b)[names(b) %in% names(RES)]
RES = b[RES, on = key]
setnames(RES, 'dif', ii)

}

}

setcolorder(RES, c('variable', pt[!grepl('Y|D', pt)], c('p', paste0(if (w) {'w_'} else {''}, tolower(pt)))))
copy(RES)
}

gbias2 = function(obs, ctrl, w = FALSE, daily_by = 'M1', wet_int_only = TRUE, fun = cor, wdt = 0.1, vars = obs[, unique(variable)], p = 0.5){

pt = attr(obs, 'terms')\$pterms
pt = pt[pt!='G1']
i = pt[10]

for (i in pt){

xo = get_value(obs, i, daily_by = daily_by, w = w)[variable %in% vars, ]
xc = get_value(ctrl, i, daily_by = daily_by, w = w)[variable %in% vars, ]
ii = if (!w) (tolower(i)) else (paste0('w_', tolower(i)))
xo = dcast(xo, formula = ... ~ variable, value.var = ii)
xc = dcast(xc, formula = ... ~ variable, value.var = ii)
setnames(xo, ncol(xo) - 1:0, c('a', 'b'))
setnames(xc, ncol(xo) - 1:0, c('a', 'b'))

if (grepl('Y|G', i)){

qxo = xo[, .(p, O = fun(a, b))]
qxc = xc[, .(p, C = fun(a, b))]
b = qxo[qxc, on = c( 'p')]
}

if (grepl('M', i)){
qxo = xo[, .(p, O = fun(a, b)), by = .(eval(parse(text = i)))]
qxc = xc[, .(p, C = fun(a, b)), by = .(eval(parse(text = i)))]
b = qxo[qxc, on = c('parse', 'p')]
}

if (grepl('D', i) & !grepl('D1', i)){
qxo = xo[, .(p, O = fun(a, b)), by = .( db = eval(parse(text = daily_by)))]
qxc = xc[, .(p, C = fun(a, b)), by = .(db = eval(parse(text = daily_by)))]
b = qxo[qxc, on = c( 'db',  'p')]
setnames(b, 'db', daily_by)
}

if (grepl('D1', i)){
qxo = xo[, .(p, O = fun(a, b)), by = .( db = eval(parse(text = daily_by)))]
qxc = xc[, .(p, C = fun(a, b)), by = .(db = eval(parse(text = daily_by)))]
b = qxo[qxc, on = c('db', 'p')]
setnames(b, 'db', daily_by)
}

b[, dif := C-O]
if ('parse' %in% names(b)) setnames(b, 'parse', i)

if ( i==pt[1] ) {
ppt = unique(c(pt[grepl('G|M', pt)], daily_by))
RES = obs[, .(p), by = ppt]
RES = b[, .(p, dif)][RES, on = c( 'p')]
setnames(RES, 'dif', if(!w) (tolower(i)) else ('w_y1') )
} else {

b[, c('O', 'C') := NULL]
key = names(b)[names(b) %in% names(RES)]
RES = b[RES, on = key]
setnames(RES, 'dif', ii)

}

}

setcolorder(RES, c( pt[!grepl('Y|D', pt)], c('p', paste0(if (w) {'w_'} else {''}, tolower(pt)))))
copy(RES)
}

melt_bias = function(b){
a1 = names(b)[names(b) %in% toupper(names(obs)) & names(b)!='variable']
a2 = names(b)[names(b) %in% tolower(names(obs)) & names(b)!='variable']

mm = melt(b, measure = a2, id = c('variable', 'p', a1))
setnames(mm, c('variable.1', 'value'), c( 'scale', 'bias'))

mm[, ts:=tscale(scale[1]), by = scale]
copy(mm)
}

melt_bias2 = function(b){
a1 = names(b)[names(b) %in% toupper(names(obs)) & names(b)!='variable']
a2 = names(b)[names(b) %in% tolower(names(obs)) & names(b)!='variable']

mm = melt(b, measure = a2, id = c('p', a1))
setnames(mm, c('variable', 'value'), c( 'scale', 'bias'))

mm[, ts:=tscale(scale[1]), by = scale]
copy(mm)
}

plot_bias = function(mb, q = 0.5, var = 'PR', nx = .1){

g = ggplot(mb[round(p, 2) == q & variable == var]) + geom_hline(yintercept = if(var=='TAS')(0)else(1), col = 'black', lty = 5)+ geom_line(aes(x = log(ts), y = bias, col = factor(M3), group = M1 ))+ scale_x_reverse(breaks = log(unique(mb\$ts)), labels = unique(mb\$scale) )+scale_y_log10()+ geom_text(aes(x = log(ts), y = bias, label = M1), data = mb[ts == min(ts) & round(p, 2) == q & variable == var], nudge_x = nx, check_overlap = TRUE)+theme_linedraw()+theme(legend.position = c(0, 1), legend.justification = c(0, 1))
structure(g)
}

plot_bias_smry = function(obs, ctrl, cor, fun = sd, keep_scale = TRUE, ...){
g1 = plot_bias(melt_bias(gbias(obs, ctrl, fun = fun, ...)), var = 'PR', q = .5) + scale_color_brewer(palette= 'Set1')+scale_y_continuous()+facet_wrap(~'PR (ORIG)')
g2 = plot_bias(melt_bias(gbias(obs, ctrl, fun = fun, ...)), var = 'TAS', q = .5) + scale_color_brewer(palette= 'Set1')+scale_y_continuous()+facet_wrap(~'TAS (ORIG)')
g3 = plot_bias(melt_bias(gbias(obs, cor, fun = fun, ...)), var = 'PR', q = .5) + scale_color_brewer(palette= 'Set1')+scale_y_continuous()+facet_wrap(~'PR (CORR)')+if(keep_scale) (ylim(range(g1\$data\$bias, na.rm = TRUE))) else (NULL)
g4 = plot_bias(melt_bias(gbias(obs, cor, fun = fun, ...)), var = 'TAS', q = .5) + scale_color_brewer(palette= 'Set1')+scale_y_continuous()+facet_wrap(~'TAS (CORR)')+if(keep_scale) (ylim(range(g2\$data\$bias, na.rm = TRUE))) else (NULL)
G = arrangeGrob(g1, g3, g2, g4, ncol = 2, nrow =2)
return(G)
}

plot_gbias = function(obs, ctrl, corr, vars = obs[, unique(variable)], fun = sd, keep_scale = TRUE, nx = 0.1, ...){

#vars = obs[, unique(variable)]
map = expand.grid('obs', c('ctrl', 'corr'), vars )
B = list()
for (i in 1:nrow(map)){
n = paste0(map[i, 'Var3'], ' (', if (map[i, 'Var2']=='ctrl')('ORIG') else ( 'CORR'), ')')
bi = gbias(obs = structure(obs[variable==map[i, 'Var3']], terms = attr(obs, 'terms')), ctrl = eval(parse(text = as.character(map[i, 'Var2']) ))[variable==map[i, 'Var3']], fun = fun, ...)
#bi = gbias(obs = obs, ctrl = eval(parse(text = as.character(map[i, 'Var2']) )), fun = fun, ...)
mbi = melt_bias(bi)
mbi[, TYPE:=if (map[i, 'Var2']=='ctrl')('ORIG') else ( 'CORR')]
B[[i]] = mbi
}
B = do.call(rbind, B)
B[, TYPE:=factor(TYPE, levels = c('ORIG', 'CORR'))]
YI = B[, .(yi= ifelse(variable[1]%in%c('TAS', 'Xd') ,0,1) ), by = .(variable, TYPE)]
G = ggplot(B[round(p, 2) == p ]) + geom_hline(aes(yintercept = yi), col = 'black', lty = 5, data = YI)+ geom_line(aes(x = log(ts), y = bias, col = factor(M3), group = M1 ))+ scale_x_reverse(breaks = log(unique(B\$ts)), labels = unique(B\$scale) )+ geom_text(aes(x = log(ts), y = bias, label = M1), data = B[ts == min(ts) & round(p, 2) == p], nudge_x = nx, check_overlap = TRUE)+theme_linedraw()+theme(legend.position = c(0, 1), legend.justification = c(0, 1))+ scale_color_brewer(palette= 'Set1') + if (keep_scale) (facet_grid(variable~TYPE, scale = 'free')) else (facet_wrap(variable~TYPE, scale = 'free', ncol = 2))

structure(G, bias = B)
}

plot_gbias2 = function(obs, ctrl, corr, fun = cor, keep_scale = TRUE, nx = 0.1, ...){

vars = obs[, unique(variable)]
cmb = t(combn(vars, 2))
map = data.table(TYP = rep(c('ctrl', 'corr'),  each = nrow(cmb)), cmb)

B = list()
for (i in 1:nrow(map)){
n = paste0(map[i, 'Var3'], ' (', if (map[i, 'Var2']=='ctrl')('ORIG') else ( 'CORR'), ')')
bi = gbias2(obs = obs, ctrl = eval(parse(text = as.character(map[i, TYP]) )), fun = fun, vars = map[i, c(V1, V2)] )
mbi = melt_bias2(bi)
mbi[, TYPE:=if (map[i, TYP]=='ctrl')('ORIG') else ( 'CORR')]
mbi[, variable:= map[i, paste(V1, V2, sep='-')]]
B[[i]] = mbi
}
B = do.call(rbind, B)
B[, TYPE:=factor(TYPE, levels = c('ORIG', 'CORR'))]
#YI = B[, .(yi= ifelse(variable[1]=='TAS',0,1) ), by = .(variable, TYPE)]
G = ggplot(B[round(p, 2) == p ]) + geom_hline(aes(yintercept = 0), col = 'black', lty = 5)+ geom_line(aes(x = log(ts), y = bias, col = factor(M3), group = M1 ))+ scale_x_reverse(breaks = log(unique(B\$ts)), labels = unique(B\$scale) )+ geom_text(aes(x = log(ts), y = bias, label = M1), data = B[ts == min(ts) & round(p, 2) == p], nudge_x = nx, check_overlap = TRUE)+theme_linedraw()+theme(legend.position = c(0, 1), legend.justification = c(0, 1))+ scale_color_brewer(palette= 'Set1') + if (keep_scale) (facet_grid(variable~TYPE, scale = 'free')) else (facet_wrap(variable~TYPE, scale = 'free', ncol = 2))

structure(G, bias = B)
}

gr4j = function(data, x, return_components = FALSE){

er=gr4j.sim(cbind(P=data\$P, E=data\$E),x1=x[1],etmult = 1, S_0 = 0.5, return_state = TRUE)
re = gr4jrouting.sim(er[, 'U'],x[2],x[3],x[4], eps=0.001, R_0 = .5, return_components = return_components )
res = cbind(er, re)
res
}

run_bil = function(dta, b = NULL, vars = 'RM'){
if (is.null(b)) {b = bil.new(type = 'd', data = dta[, .(DTM, P = PR, T = TAS)]) } else {
bil.set.values(b, dta[, .(DTM, P = PR, T = TAS)])
}
bil.pet(b)
d = data.frame(data.frame(bil.run(b))[, vars])
names(d) = vars
d
}

run_gr4j = function(dta, g = NULL, vars = c('U', 'S', 'ET', 'Xr', 'Xd', 'R', 'Ps', 'Perc', 'Es')){
b = bil.new(type = 'd', data = dta[, .(DTM, P = PR, T = TAS)])
bil.pet(b)
d = bil.get.data(b)
#rm1 = gr4j(data.frame(d[, .(P, E = PET)]), .g\$optim\$optim\$bestmem, return_components = FALSE)
rm = gr4j(data.frame(d[, .(P, E = PET)]), .g\$optim\$optim\$bestmem, return_components = TRUE)
rm[1, c('Xr', 'Xd')] = 0
X1 = .g\$optim\$optim\$bestmem[1]
rm = data.table(DTM = d\$DTM, rm)
rm[, P := d\$P]
rm[, E := d\$PET]
rm[, Pn := max(0, P-E), by = DTM]
rm[, En := max(0, E-P), by = DTM]
rm[, Ps := (X1 * ( 1 - (S / X1) ^2 ) * tanh(Pn/X1)) / (1 + S/X1 * tanh(Pn/X1))]
rm[, Es := (S * (2 - S/X1) * tanh(En/X1) )/(1 + (1 - S/X1) * tanh(En/X1))]
rm[, Perc := S * (1 - (1 + (4/9 * S/X1)^4 ) ^ (1/4) )]
# #rm = data.frame(rm1, rm2)
# colnames(rm) =
# #data.frame(rm)[, vars]
copy(rm)
}
```
hanel/muscabicor documentation built on May 14, 2017, 10:25 p.m.