svyafc <-
  function(formula,
           design,
           k ,
           g ,
           cutoffs ,
           dimw = NULL,
           na.rm = FALSE,
           ...) {
    if (k <= 0 |
        k > 1)
      stop("This functions is only defined for k in (0,1].")
    if (g < 0)
      stop("This function is undefined for g < 0.")
    if (!is.list(cutoffs))
      stop("The parameter 'cutoffs' has to be a list.")
    
    ach.matrix <-
      model.frame(formula, design$variables, na.action = na.pass)[, ]
    
    if (!is.null(dimw)) {
      if (any(is.na(dimw))) {
        stop("Invalid value in dimension weights vector.")
      }
      if (sum(dimw) > 1) {
        stop("The sum of dimension weigths have to be equal to one.")
      }
      if (any(dimw > 1 |
              dimw < 0)) {
        stop("Dim. weights have to be within interval [0,1].")
      }
      if (length(dimw) != ncol(ach.matrix)) {
        stop("Dimension weights' length differs from number of dimensions in formula")
      }
    }
    
    var.class <- lapply(ach.matrix, function(x)
      class(x)[1])
    var.class <- matrix(
      var.class,
      nrow = 1,
      ncol = ncol(ach.matrix),
      dimnames = list(c("var.class"), colnames(ach.matrix))
    )
    
    if (any(!(var.class %in% c("numeric", "integer", "ordered")))) {
      stop(
        "This function is only applicable to variables of types 'numeric' or 'ordered factor'."
      )
    }
    if (any((var.class == "integer"))) {
      stop(
        "At least one of the variables is an integer.\nCoerce your column to numeric with as.numeric if you are sure it's what you want."
      )
    }
    
    w <- 1 / design$prob
    
    if (any(ach.matrix[w != 0, var.class == "numeric"] < 0, na.rm = TRUE))
      stop(
        "The Alkire-Foster multidimensional poverty class is defined for non-negative numeric variables only."
      )
    
    if (na.rm) {
      nas <- apply(ach.matrix, 1, function(x)
        any(is.na(x)))
      design <- design[nas == 0,]
      w <- 1 / design$prob
    }
    
    ach.matrix <-
      model.frame(formula, design$variables, na.action = na.pass)[, ]
    ach.matrix <- ach.matrix [w > 0,]
    w <- w [w > 0]
    
    if (any(is.na(ach.matrix))) {
      if (is.null(dimw)) {
        dimw = rep(1 / ncol(var.class), length(var.class))
      }
      
      rval <- as.numeric(NA)
      variance <- as.numeric(NA)
      class(rval) <- c("cvystat" , "svystat")
      attr(rval, "var") <- variance
      names(rval)[1] <- attr(rval, "statistic") <- "alkire-foster"
      dimtable <-
        as.data.frame(matrix(
          c(strsplit(as.character(formula)[[2]] , ' \\+ ')[[1]], dimw),
          nrow = ncol(var.class),
          ncol = 2,
          dimnames = list(paste("dimension", 1:ncol(var.class)), c("variables", "weight"))
        ),
        stringsAsFactors = FALSE)
      dimtable[, 2] <- as.numeric(dimtable[, 2])
      attr(rval, "dimensions") <- dimtable
      attr(rval, "parameters") <-
        matrix(
          c(g, k),
          nrow = 1,
          ncol = 2,
          dimnames = list("parameters", c("g=", "k="))
        )
      
      return(rval)
      
    }
    
    # Deprivation Matrix
    dep.matrix <-
      sapply(
        seq_along(var.class),
        FUN = function(i) {
          1 * (cutoffs[[i]] > ach.matrix[, i])
        }
      )
    colnames(dep.matrix) <- colnames(var.class)
    
    # Unweighted count of deprivations:
    # depr.count <- rowSums( dep.matrix )
    
    # deprivation k cut
    if (is.null(dimw)) {
      dimw = rep(1 / ncol(dep.matrix), ncol(dep.matrix))
    }
    
    # Weighted sum of deprivations:
    depr.sums <- rowSums(sweep(dep.matrix, MARGIN = 2 , dimw, `*`))
    
    # k multidimensional cutoff:
    multi.cut <- depr.sums * (depr.sums >= k)
    rm(dep.matrix)
    
    # Censored Deprivation Matrix
    cen.dep.matrix <-
      sapply(
        seq_along(cutoffs) ,
        FUN = function(x) {
          if (var.class[x] == "numeric") {
            1 * (cutoffs[[x]] > ach.matrix[, x]) * ((cutoffs[[x]] - ach.matrix[, x]) / cutoffs[[x]]) ^
              g
          } else {
            1 * (cutoffs[[x]] > ach.matrix[, x])
          }
        }
      )
    colnames(cen.dep.matrix) <- colnames(var.class)
    cen.dep.matrix[multi.cut == 0,] <- 0
    
    # Sum of censored deprivations:
    cen.depr.sums <-
      rowSums(sweep(cen.dep.matrix, MARGIN = 2 , dimw, `*`))
    rm(cen.dep.matrix, ach.matrix)
    
    
    w <- 1 / design$prob
    w[w > 0] <- cen.depr.sums
    cen.depr.sums <- w
    rm(w)
    
    if (g == 0) {
      w <- 1 / design$prob
      w[w > 0] <- (depr.sums >= k)
      h_i <- w
      rm(w)
      
      h_est <- survey::svymean(h_i, design)
      
      w <- 1 / design$prob
      w[w > 0] <- multi.cut
      multi.cut <- w
      rm(w)
      
      a_est <- survey::svyratio(multi.cut, h_i, design)
      
      
    }
    
    rval <- survey::svymean(cen.depr.sums , design)
    names(rval)[1] <- attr(rval, "statistic") <- "alkire-foster"
    dimtable <-
      as.data.frame(matrix(
        c(strsplit(as.character(formula)[[2]] , ' \\+ ')[[1]], dimw),
        nrow = ncol(var.class),
        ncol = 2,
        dimnames = list(paste("dimension", 1:ncol(var.class)), c("variables", "weight"))
      ),
      stringsAsFactors = FALSE)
    dimtable[, 2] <- as.numeric(dimtable[, 2])
    attr(rval, "dimensions") <- dimtable
    attr(rval, "parameters") <-
      matrix(
        c(g, k),
        nrow = 1,
        ncol = 2,
        dimnames = list("parameters", c("g=", "k="))
      )
    if (g == 0) {
      attr(rval, "extra") <-
        matrix(
          c(h_est[1], a_est[[1]], attr(h_est, "var")[1] ^ .5, a_est[[2]] ^ .5),
          nrow = 2,
          ncol = 2,
          dimnames = list(c("H", "A"), c("coef", "SE"))
        )
    }
    class(rval) <- c("cvystat" , "svystat")
    
    return(rval)
    
  }