bael_calib <- function(formula = waga ~ plec:klm2:grwiek1 | idm + lp | msc,
                       data,
                       totals,
                       method = 'linear',
                       bounds = c(0,10)) {

  fm1 <- as.Formula(formula)

  vars_weight <- attr(fm1, 'lhs')[[1]]
  vars_calib <- attr(fm1, 'rhs')[[1]]
  vars_ids <- attr(fm1, 'rhs')[[2]]
  vars_split <- as.character(attr(fm1, 'rhs')[[3]])


  vars_calib <- unlist(strsplit(paste(vars_calib), '\\+|\\*|:'))
  vars_calib <- vars_calib[nchar(vars_calib) > 0]

  vars_ids <- unlist(strsplit(paste(vars_ids), '\\+|\\*|:'))
  vars_ids <- vars_ids[nchar(vars_ids) > 0]

  ## tworzenie zmiennych 0-1 -- na podstawie
  setkeyv(data, c(vars_split)) ## dla przyspieszenia

  data[, aux_x := Reduce(function(...) paste0(...), .SD), .SDcols = vars_calib]
  data[, ids := Reduce(function(...) paste0(...), .SD), .SDcols = vars_ids]

  ## dane do kalibracji
  aux <- data[, list(aux = list(t(fac2sparse(aux_x))),
                     ids = list(ids)),
              by = c(vars_split)]

  d <- data[, list(d = list(as.numeric(unlist(.SD)))),
            by = c(vars_split),
            .SDcols = as.character(vars_weight)]

  totals <- totals[, list(pop = list(pop)),
                   by = c(vars_split)]

  do_kalibracji <- totals[d[aux,
                          on = vars_split],
                          on = vars_split]
  ## test czy długości totali i wymiary X są równe
  do_kalibracji[,flag:=length(pop[[1]])==ncol(aux[[1]]),
                by = c(vars_split)]

  method_no <- length(method)

  if (method_no == 1) {
     po_kalib <- do_kalibracji[flag == TRUE,
                            list(g = list(calibWeights2(d = d[[1]],
                                                        X = aux[[1]],
                                                        totals = pop[[1]],
                                                        method = method,
                                                        bounds = bounds)),
                                 ids = list(ids[[1]])),
    by = c(vars_split)]

    result <- po_kalib[,list(g = unlist(g),
                 ids = unlist(ids)),
                 by = c(vars_split)]
  } else {
    po_kalib <- list()
    for (m in method) {
        po_kalib[[m]] <- do_kalibracji[flag == TRUE,
                            list(g = list(calibWeights2(d = d[[1]],
                                                        X = aux[[1]],
                                                        totals = pop[[1]],
                                                        method = m,
                                                        bounds = bounds)),
                                 ids = list(ids[[1]])),
                            by = c(vars_split)]

        po_kalib_ <- lapply(po_kalib, function(x) {
         x[,list(g = unlist(g),
                 ids = unlist(ids)),
           by = c(vars_split)]
       })
      result <- rbindlist(po_kalib_,idcol = 'method')
    }
  }
  return(result)

}