R/decomps.R

Defines functions run_gr4j run_bil gr4j plot_gbias2 plot_gbias plot_bias_smry plot_bias melt_bias2 melt_bias gbias2 gbias bias get_value give_map decomp rev_dif dif agg prob

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.