R/family.basics.R

In VGAM: Vector Generalized Linear and Additive Models

# These functions are
# Copyright (C) 1998-2024 T.W. Yee, University of Auckland.
# All rights reserved.













getarg <- function(a) {



  b <- unlist(strsplit(a, "(", fixed = TRUE))

  if (length(b) == 2) {  # Usual case
    d <- unlist(strsplit(b[2], ")", fixed = TRUE))
    d[1]
  } else if (length(b) == 1) {  # identitylink
    b
  } else if (length(b) == 3) {  # nbcanlink
    d <- unlist(strsplit(b[2], ",", fixed = TRUE))
    d[1]
  } else {
    a
  }
}  # getarg








subsetcol <-
Select <-
  function(
           data = list(),
           prefix = "y",

           lhs  = NULL,
           rhs  = NULL,  # Can be "0" to suppress an intercept, else "".
           rhs2 = NULL,  # Can be "0" to suppress an intercept, else "".
           rhs3 = NULL,  # Can be "0" to suppress an intercept, else "".

           as.character = FALSE,
           as.formula.arg = FALSE,
           tilde = TRUE,
           exclude = NULL,
           sort.arg = TRUE) {





  if (is.character(exclude))
    if (any(nchar(prefix) == 0))
      stop("bad input for argument 'exclude'")
  if (!is.logical(sort.arg) ||
      length(sort.arg) != 1)
    stop("bad input for argument 'sort.arg'")


  col.names <- colnames(data)
  if (is.logical(prefix)) {
    index <- if (prefix) seq_along(col.names) else
             stop("cannot have 'prefix = FALSE'")
  } else {
    index <- NULL
    for (ii in seq_along(prefix)) {
      small.col.names <- substr(col.names, 1, nchar(prefix[ii]))

      index <- c(index, grep(prefix[ii], small.col.names))
    }
  }





  temp.col.names <- col.names[index]


  if (length(exclude)) {
    exclude.index <- NULL
    for (ii in seq_along(exclude)) {
      exclude.index <- c(exclude.index,
                  (seq_along(col.names))[exclude[ii] == col.names])
    }
    exclude.index <- unique(sort(exclude.index))
    index <- setdiff(index, exclude.index)
    temp.col.names <- col.names[index]
  }




  if (sort.arg) {
    ooo <- order(temp.col.names)
    index <- index[ooo]
    temp.col.names <- temp.col.names[ooo]
  }



 ltcn.positive <- (length(temp.col.names) > 0)



  if (as.formula.arg) {
    form.string <- paste0(ifelse(length(lhs), lhs, ""),
                          ifelse(tilde, " ~ ", ""),
                          if (ltcn.positive)
                            paste(temp.col.names, collapse = " + ") else
                            "",
                      ifelse(ltcn.positive && length(rhs ), " + ", ""),
                          ifelse(length(rhs ), rhs, ""),
                          ifelse(length(rhs2), paste(" +", rhs2), ""),
                          ifelse(length(rhs3), paste(" +", rhs3), ""))

    if (as.character) {
      form.string
    } else {
      as.formula(form.string)
    }
  } else {
    if (as.character) {
      paste0("cbind(",
             paste(temp.col.names, collapse = ", "),
             ")")
    } else {
      ans <- if (is.matrix(data)) data[, index] else
             if (is.list(data)) data[index] else
             stop("argument 'data' is neither a list or a matrix")
      if (length(ans)) {
        as.matrix(ans)
      } else {
        NULL
      }
    }
  }
}  # subsetcol, Select








if (FALSE)
subsetc <-
  function(x, select,
           prefix = NULL,
           subset = TRUE, drop = FALSE,
           exclude = NULL,
           sort.arg = !is.null(prefix),
           as.character = FALSE) {

  if (!is.null(prefix)) {
    if (!missing(select))
      warning("overwriting argument 'select' by something ",
              "using 'prefix'")
    select <- grepl(paste0("^", prefix), colnames(x))
  }

  if (missing(select)) {
    vars <- TRUE
  } else {
    nl <- as.list(seq_along(x))  # as.list(1L:ncol(x))
    names(nl) <- names(x)  # colnames(x)
    vars <- eval(substitute(select), nl, parent.frame())
  }

  ans <- x[subset & !is.na(subset), vars, drop = drop]
  if (sort.arg) {
    cna <- colnames(ans)
    ooo <- order(cna)
    ans <- ans[, ooo, drop = drop]
  }

  if (!is.null(exclude)) {
    cna <- colnames(ans)
    ooo <- match(exclude, cna)
    ans <- ans[, -ooo, drop = drop]
  }

  if (as.character) {
    cna <- colnames(ans)
    paste0("cbind(", paste0(cna, collapse = ", "), ")")
  } else {
    ans
  }
}







 grid.search <- function(vov, objfun, y, x, w, extraargs = NULL,
                         maximize = TRUE, abs.arg = FALSE,
                         ret.objfun = FALSE, ...) {
  if (!is.vector(vov))
    stop("argument 'vov' must be a vector")
  objvals <- vov
  for (ii in seq_along(vov))
    objvals[ii] <- objfun(vov[ii], y = y, x = x, w = w,
                          extraargs = extraargs, ...)
  try.this <- if (abs.arg) {
               if (maximize) vov[abs(objvals) == max(abs(objvals))] else
               vov[abs(objvals) == min(abs(objvals))]
             } else {
               if (maximize) vov[objvals == max(objvals)] else
               vov[objvals == min(objvals)]
             }
  if (!length(try.this))
    stop("something has gone wrong!")
  ans <- if (length(try.this) == 1)
    try.this else sample(try.this, size = 1)


  myvec <- objvals[ans == vov]  # Could be a vector
  if (ret.objfun) c(Value = ans, ObjFun = myvec[1]) else ans
}  # grid.search






 grid.search2 <-
    function(vov1, vov2, objfun, y, x, w, extraargs = NULL,
             maximize = TRUE,  # abs.arg = FALSE,
             ret.objfun = FALSE, ...) {
  if (!is.vector(vov1))
    stop("argument 'vov1' must be a vector")
  if (!is.vector(vov2))
    stop("argument 'vov2' must be a vector")

  allmat1 <- expand.grid(vov1 = as.vector(vov1),
                         vov2 = as.vector(vov2))
  objvals <- numeric(nrow(allmat1))
  for (ii in seq_along(objvals))
    objvals[ii] <- objfun(allmat1[ii, "vov1"],
                          allmat1[ii, "vov2"],
                          y = y, x = x, w = w,
                          extraargs = extraargs, ...)

  ind5 <- if (maximize) which.max(objvals) else which.min(objvals)


  c(Value1 = allmat1[ind5, "vov1"],
    Value2 = allmat1[ind5, "vov2"],
    ObjFun = if (ret.objfun) objvals[ind5] else NULL)
}  # grid.search2




 grid.search3 <-
    function(vov1, vov2, vov3, objfun, y, x, w, extraargs = NULL,
             maximize = TRUE,  # abs.arg = FALSE,
             ret.objfun = FALSE, ...) {
  if (!is.vector(vov1))
    stop("argument 'vov1' must be a vector")
  if (!is.vector(vov2))
    stop("argument 'vov2' must be a vector")
  if (!is.vector(vov3))
    stop("argument 'vov3' must be a vector")

  allmat1 <- expand.grid(vov1 = as.vector(vov1),
                         vov2 = as.vector(vov2),
                         vov3 = as.vector(vov3))
  objvals <- numeric(nrow(allmat1))
  for (ii in seq_along(objvals))
    objvals[ii] <- objfun(allmat1[ii, "vov1"],
                          allmat1[ii, "vov2"],
                          allmat1[ii, "vov3"],
                          y = y, x = x, w = w,
                          extraargs = extraargs, ...)

  ind5 <- if (maximize) which.max(objvals) else which.min(objvals)


  c(Value1 = allmat1[ind5, "vov1"],
    Value2 = allmat1[ind5, "vov2"],
    Value3 = allmat1[ind5, "vov3"],
    ObjFun = if (ret.objfun) objvals[ind5] else NULL)
}  # grid.search3






 grid.search4 <-
    function(vov1, vov2, vov3, vov4,
             objfun, y, x, w, extraargs = NULL,
             maximize = TRUE,  # abs.arg = FALSE,
             ret.objfun = FALSE, ...) {
  if (!is.vector(vov1))
    stop("argument 'vov1' must be a vector")
  if (!is.vector(vov2))
    stop("argument 'vov2' must be a vector")
  if (!is.vector(vov3))
    stop("argument 'vov3' must be a vector")
  if (!is.vector(vov4))
    stop("argument 'vov4' must be a vector")

  allmat1 <- expand.grid(vov1 = as.vector(vov1),
                         vov2 = as.vector(vov2),
                         vov3 = as.vector(vov3),
                         vov4 = as.vector(vov4))
  objvals <- numeric(nrow(allmat1))
  for (ii in seq_along(objvals))
    objvals[ii] <- objfun(allmat1[ii, "vov1"],
                          allmat1[ii, "vov2"],
                          allmat1[ii, "vov3"],
                          allmat1[ii, "vov4"],
                          y = y, x = x, w = w,
                          extraargs = extraargs, ...)

  ind5 <- if (maximize) which.max(objvals) else which.min(objvals)


  c(Value1 = allmat1[ind5, "vov1"],
    Value2 = allmat1[ind5, "vov2"],
    Value3 = allmat1[ind5, "vov3"],
    Value4 = allmat1[ind5, "vov4"],
    ObjFun = if (ret.objfun) objvals[ind5] else NULL)
}  # grid.search4






 getind <- function(constraints, M, ncolx) {



  if (!length(constraints)) {

    constraints <- vector("list", ncolx)
    for (ii in 1:ncolx)
      constraints[[ii]] <- diag(M)
  }

  ans <- vector("list", M+1)
  names(ans) <- c(param.names("eta", M), "ncolX.vlm")

  temp2 <- matrix(unlist(constraints), nrow = M)
  for (kk in 1:M) {
    ansx <- NULL
    for (ii in seq_along(constraints)) {
      temp <- constraints[[ii]]
      isfox <- any(temp[kk, ] != 0)
      if (isfox) {
        ansx <- c(ansx, ii)
      }
    }
    ans[[kk]] <- list(xindex = ansx,
                  X.vlmindex = (1:ncol(temp2))[temp2[kk,] != 0])
  }
  ans[[M+1]] <- ncol(temp2)

  ans
}  # genind




 cm.VGAM <-
  function(cm, x, bool, constraints,
           apply.int = FALSE,
           cm.default = diag(nrow(cm)),  # 20121226
           cm.intercept.default = diag(nrow(cm))  # 20121226
          ) {




  if (is.null(bool))
    return(NULL)

  if (!is.matrix(cm))
    stop("argument 'cm' is not a matrix")
  M <- nrow(cm)
  asgn <- attr(x, "assign")
  if (is.null(asgn))
    stop("the 'assign' attribute is missing from 'x'; this ",
         "may be due to some missing values")  # 20100306
  nasgn <- names(asgn)
  ninasgn <- nasgn[nasgn != "(Intercept)"]

  if (!length(constraints)) {
    constraints <- vector("list", length(nasgn))
    for (ii in seq_along(nasgn)) {
      constraints[[ii]] <- cm.default  # diag(M)
    }
    names(constraints) <- nasgn


    if (any(nasgn == "(Intercept)"))
      constraints[["(Intercept)"]] <- cm.intercept.default
  }

  if (!is.list(constraints))
    stop("argument 'constraints' must be a list")

  if (length(constraints) != length(nasgn) ||
      any(sort(names(constraints)) != sort(nasgn))) {
    cat("\nnames(constraints)\n")
    print(names(constraints) )
    cat("\nnames(attr(x, 'assign'))\n")
    print( nasgn )
    stop("The above do not match; 'constraints' is half-pie")
  }




  if (is.logical(bool)) {
    if (bool) {
      if (any(nasgn == "(Intercept)") && apply.int)
        constraints[["(Intercept)"]] <- cm


      if (length(ninasgn))
        for (ii in ninasgn)
          constraints[[ii]] <- cm
    } else {
      return(constraints)
    }
  } else {
    tbool <- terms(bool)
    if (attr(tbool, "response")) {
      ii <- attr(tbool, "factors")
      default <- dimnames(ii)[[1]]
      default <- default[1]
      default <- if (is.null(default[1])) {
        t.or.f <- attr(tbool, "variables")

        t.or.f <- as.character( t.or.f )
        if (t.or.f[1] == "list" && length(t.or.f) == 2 &&
           (t.or.f[2] == "TRUE" || t.or.f[2] == "FALSE")) {
          t.or.f <- as.character( t.or.f[2] )
          parse(text = t.or.f)[[1]]
        } else {
          stop("something gone awry")
        }
      } else {
        parse(text = default[1])[[1]]  # Original
      }
      default <- as.logical(eval(default))
    } else {
      default <- TRUE
    }
    tl <- attr(tbool, "term.labels")
    if (attr(tbool, "intercept"))
      tl <- c("(Intercept)", tl)

    for (ii in nasgn) {
      if ( default &&  any(tl == ii))
        constraints[[ii]] <- cm
      if (!default && !any(tl == ii))
        constraints[[ii]] <- cm
    }
  }

  constraints
}  # cm.VGAM






cm.nointercept.VGAM <- function(constraints, x, nointercept, M) {

  asgn <- attr(x, "assign")
  nasgn <- names(asgn)
  if (is.null(constraints)) {
    constraints <- vector("list", length(nasgn))  # list()
    names(constraints) <- nasgn
  }
  if (!is.list(constraints))
    stop("'constraints' must be a list")

  for (ii in seq_along(asgn))
    constraints[[nasgn[ii]]] <- if (is.null(constraints[[nasgn[ii]]]))
      diag(M) else eval(constraints[[nasgn[ii]]])

  if (is.null(nointercept))
    return(constraints)
  if (!is.numeric(nointercept))
    stop("'nointercept' must be numeric")

  nointercept <- unique(sort(nointercept))
  if (length(nointercept) == 0 || length(nointercept) >= M)
    stop("too few or too many values")

  if (any(nointercept < 1 | nointercept > M))
    stop("'nointercept' out of range")
  if (nasgn[1] != "(Intercept)" || M == 1)
    stop("Need an (Intercept) constraint matrix with M>1")
  if (!identical(constraints[["(Intercept)"]], diag(M)))
    warning("Constraint matrix of (Intercept) not diagonal")

  temp <- constraints[["(Intercept)"]]
  temp <- temp[, -nointercept, drop = FALSE]
  constraints[["(Intercept)"]] <- temp
  constraints
}  # cm.nointercept.VGAM






 cm.zero.VGAM <- function(constraints, x, zero = NULL, M = 1,
                          predictors.names, M1 = 1) {



  if (is.character(predictors.names)) {
    for (ii in 1:length(predictors.names))
      predictors.names[ii] <- getarg(predictors.names[ii])
  }




  dotzero <- zero  # Transition

  if (length(dotzero) == 1 && (dotzero == "" || is.na(dotzero)))
    dotzero <- NULL

  if (is.character(dotzero)) {




    which.numeric.all <- NULL
    for (ii in seq_along(dotzero)) {
      which.ones <-
          grep(dotzero[ii], predictors.names, fixed = TRUE)
      if (length(which.ones)) {
        which.numeric.all <- c(which.numeric.all, which.ones)
      } else {
        warning("some values of argument 'zero' are unmatched. ",
                "Ignoring them")
      }
    }  # for ii
    which.numeric <- unique(sort(which.numeric.all))

    if (!length(which.numeric)) {
      warning("No values of argument 'zero' were matched.")
      which.numeric <- NULL
    } else if (length(which.numeric.all) > length(which.numeric)) {
      warning("There were redundant values of argument 'zero'.")
    }

    dotzero <- which.numeric
  }  # if is.character(dotzero)



  posdotzero <-  dotzero[dotzero > 0]
  negdotzero <-  dotzero[dotzero < 0]


  zneg.index <- if (length(negdotzero)) {

    if (!is.Numeric(-negdotzero, positive = TRUE,
                    integer.valued = TRUE) ||
        max(-negdotzero) > M1)
        stop("bad input for argument 'zero'")

    bigUniqInt <- 1080
    zneg.index <- rep(0:bigUniqInt, rep(length(negdotzero),
                      1 + bigUniqInt)) * M1 + abs(negdotzero)
    sort(intersect(zneg.index, 1:M))
  } else {
    NULL
  }

  zpos.index <- if (length(posdotzero)) posdotzero else NULL
  z.Index <- if (!length(dotzero)) NULL else
               unique(sort(c(zneg.index, zpos.index)))


  zero <- z.Index  # Transition



  asgn <- attr(x, "assign")
  nasgn <- names(asgn)
  if (is.null(constraints)) {
    constraints <- vector("list", length(nasgn))  # list()
    names(constraints) <- nasgn
  }
  if (!is.list(constraints))
    stop("'constraints' must be a list")

  for (ii in seq_along(asgn))
    constraints[[nasgn[ii]]] <- if (is.null(constraints[[nasgn[ii]]]))
      diag(M) else eval(constraints[[nasgn[ii]]])

  if (is.null(zero))
    return(constraints)

  if (any(zero < 1 | zero > M))
    stop("argument 'zero' out of range; should have values between ",
         "1 and ", M, " inclusive")
  if (nasgn[1] != "(Intercept)")
    stop("cannot fit an intercept to a no-intercept model")

  if (2 <= length(constraints))
    for (ii in 2:length(constraints)) {
      Hmatk <- constraints[[nasgn[ii]]]
      Hmatk[zero, ] <- 0
      index <- NULL
      for (kk in 1:ncol(Hmatk))
        if (all(Hmatk[, kk] == 0))
          index <- c(index, kk)
      if (length(index) == ncol(Hmatk))
        stop("constraint matrix has no columns!")
      if (!is.null(index))
        Hmatk <- Hmatk[, -index, drop = FALSE]
      constraints[[nasgn[ii]]] <- Hmatk
    }  # for ii

  constraints
}  # cm.zero.VGAM






 process.constraints <-
  function(constraints, x, M,
           by.col = TRUE, specialCM = NULL,
           Check.cm.rank = TRUE) {







    asgn <- attr(x, "assign")
    nasgn <- names(asgn)

  if (is.null(constraints)) {
    constraints <- vector("list", length(nasgn))
    for (ii in seq_along(nasgn))
      constraints[[ii]] <- diag(M)
    names(constraints) <- nasgn
  }

  if (is.matrix(constraints))
    constraints <- list(constraints)

  if (!is.list(constraints))
    stop("'constraints' must be a list")

  lenconstraints <- length(constraints)
  if (lenconstraints > 0)
    for (ii in 1:lenconstraints) {
      list.elt <- constraints[[ii]]

      if (is.function(list.elt)) {
        list.elt <- list.elt()
      }

      constraints[[ii]] <- eval(list.elt)
      if (!is.null  (constraints[[ii]]) &&
          !is.matrix(constraints[[ii]]))
        stop("'constraints[[", ii, "]]' is not a matrix")
    }  # for ii

  if (is.null(names(constraints)))
    names(constraints) <- rep_len(nasgn, lenconstraints)

  tmp3 <- vector("list", length(nasgn))
  names(tmp3) <- nasgn
  for (ii in seq_along(nasgn))
    tmp3[[nasgn[ii]]] <-
      if (is.null(constraints[[nasgn[ii]]])) diag(M) else
             eval(constraints[[nasgn[ii]]])

  for (ii in seq_along(asgn)) {
    if (!is.matrix(tmp3[[ii]])) {
      stop("component ", ii, "is not a constraint matrix")
    }
    if (ncol(tmp3[[ii]]) > M)
      stop("constraint matrix has too many columns")
  }

  if (!by.col)
    return(tmp3)

  constraints <- tmp3
  Hlist <- vector("list", ncol(x))
  for (ii in seq_along(asgn)) {
    cols <- asgn[[ii]]
    ictr <- 0
    for (jay in cols) {
      ictr <- ictr + 1
      cm <- if (is.list(specialCM) &&
                any(nasgn[ii] == names(specialCM))) {
              slist <- specialCM[[(nasgn[ii])]]
              slist[[ictr]]
            } else {
              constraints[[ii]]
            }
      Hlist[[jay]] <- cm
    }
  }
  names(Hlist) <- colnames(x)  # dimnames(x)[[2]]



  if (Check.cm.rank) {
    all.svd.d <- function(x) svd(x)$d
    mylist <- lapply(Hlist, all.svd.d)

    if (max(unlist(lapply(mylist, length))) > M)
      stop("some constraint matrices have more than ", M,
           "columns")

    MyVector <- unlist(mylist)
    if (min(MyVector) < 1.0e-10)
      stop("some constraint matrices are not of ",
           "full column-rank: ",
           paste(names(MyVector)[MyVector < 1.0e-10], collapse = ", "))
  }

  Hlist
}  # process.constraints






 trivial.constraints <- function(Hlist, target = diag(M)) {


  if (is.null(Hlist))
    return(1)

  if (is.matrix(Hlist))
    Hlist <- list(Hlist)
  M <- dim(Hlist[[1]])[1]

  if (!is.matrix(target))
    stop("target is not a matrix")
  dimtar <- dim(target)

  trivc <- rep_len(1, length(Hlist))
  names(trivc) <- names(Hlist)
  for (ii in seq_along(Hlist)) {
    d <- dim(Hlist[[ii]])
    if (d[1] != dimtar[1]) trivc[ii] <- 0
    if (d[2] != dimtar[2]) trivc[ii] <- 0
    if (d[1] != M)         trivc[ii] <- 0
    if (length(Hlist[[ii]]) != length(target))
      trivc[ii] <- 0
    if (trivc[ii] == 0) next
    if (!all(c(Hlist[[ii]]) == c(target)))
      trivc[ii] <- 0
    if (trivc[ii] == 0) next
  }
  trivc
}  # trivial.constraints






 add.constraints <- function(constraints, new.constraints,
                             overwrite = FALSE, check = FALSE) {

  empty.list <- function(l)
    (is.null(l) || (is.list(l) && length(l) == 0))

  if (empty.list(constraints))
    if (is.list(new.constraints))
      return(new.constraints) else
      return(list())  # Both NULL probably

  constraints <- as.list(constraints)
  new.constraints <- as.list(new.constraints)
  nc <- names(constraints)         # May be NULL
  nn <- names(new.constraints)     # May be NULL

  if (is.null(nc) || is.null(nn))
    stop("lists must have names")
  if (any(nc == "") || any(nn == ""))
    stop("lists must have names")

  if (!empty.list(constraints) && !empty.list(new.constraints)) {
    for (ii in nn) {
      if (any(ii == nc)) {
        if (check &&
        (!(all(dim(constraints[[ii]]) == dim(new.constraints[[ii]])) &&
           all(    constraints[[ii]]  ==     new.constraints[[ii]]))))
          stop("apparent contradiction in the specification ",
               "of the constraints")
        if (overwrite)
          constraints[[ii]] <- new.constraints[[ii]]
      } else
        constraints[[ii]] <- new.constraints[[ii]]
    }
  } else {
    if (!empty.list(constraints))
      return(as.list(constraints)) else
      return(as.list(new.constraints))
  }

  constraints
}  # add.constraints











 iam <- function(j, k, M,  # hbw = M,
                 both = FALSE, diag = TRUE) {



  jay <- j
  kay <- k

  if (M == 1)
    if (!diag) stop("cannot handle M == 1 and diag = FALSE")

  if (M == 1) {
    if (both) return(list(row.index = 1, col.index = 1)) else return(1)
  }
     
  upper <- if (diag) M else M - 1
  i2 <- as.list(upper:1)
  i2 <- lapply(i2, seq)
  i2 <- unlist(i2)

  i1 <- matrix(1:M, M, M)
  i1 <- if (diag) c(i1[row(i1) >= col(i1)]) else
                  c(i1[row(i1) >  col(i1)])

  if (both) {
    list(row.index = i2, col.index = i1)
  } else {
    if (any(jay > M) || any(kay > M) || any(jay < 1) || any(kay < 1))
      stop("range error in argument 'j' or 'k'")

    if (length(jay) < length(kay)) jay <- rep_len(jay, length(kay))
    if (length(kay) < length(jay)) kay <- rep_len(kay, length(jay))

    if (all(c(length(jay), length(kay)) == 1)) {
      both <- (i1 == jay & i2 == kay) | (i1 == kay & i2 == jay)
      return((seq_along(i2))[both])
    } else {
    vector.big <- 10000 * pmin(i1, i2) + pmax(i1, i2)
    vector.sml <- 10000 * pmin(jay, kay) + pmax(jay, kay)
      match.sml <- match(vector.sml, vector.big)
      match.sml
    }
  }
}  # iam







if (FALSE)
 miam <- function(j, k, M1,  # M1 used to be called M
                  NOS = 1,  # This argument is new
                  both = FALSE, diag = TRUE) {

  if (NOS == 1)
    return(iam(j = j, k = k, M = M1, both = both, diag = diag))

  M <- M1 * NOS


  jay <- j
  kay <- k

  if (M == 1)
    if (!diag) stop("cannot handle this")

  if (M == 1)
    if (both) return(list(row.index = 1, col.index = 1)) else
    return(1)

  upper <- if (diag) M else M - 1
  i2 <- as.list(upper:1)
  i2 <- lapply(i2, seq)
  i2 <- unlist(i2)

  i1 <- matrix(1:M, M, M)
  i1 <- if (diag) c(i1[row(i1) >= col(i1)]) else
                  c(i1[row(i1) >  col(i1)])

  if (both) {
    list(row.index = i2, col.index = i1)
  } else {
    if (jay > M || kay > M || jay < 1 || kay < 1)
      stop("range error in j or k")
    both <- (i1 == jay & i2 == kay) |
            (i1 == kay & i2 == jay)
    (seq_along(i2))[both]
  }
}  # miam






 dimm <- function(M, hbw = M) {

  if (!is.numeric(hbw))
    hbw <- M

  if (hbw > M || hbw < 1)
    stop("range error in argument 'hbw'")
  hbw * (2 * M - hbw +1) / 2
}  # dimm









 m2a <- function(m, M, upper = FALSE, allow.vector = FALSE) {


  if (!is.numeric(m))
      stop("argument 'm' is not numeric")

  if (!is.matrix(m))
    m <- cbind(m)
  n <- nrow(m)
  dimm <- ncol(m)
  index <- iam(NA, NA, M = M, both = TRUE, diag = TRUE)
  if (dimm > length(index$row.index))
    stop("bad value for 'M'; it is too small")
  if (dimm < M) {
    stop("bad value for 'M'; it is too big")
  }

  fred <- .C("m2accc", as.double(t(m)), ans=double(M*M*n),
      as.integer(dimm),
      as.integer(index$row-1),
      as.integer(index$col-1),
      as.integer(n),  as.integer(M),
      as.integer(as.numeric(upper)), NAOK = TRUE)
  dim(fred$ans) <- c(M, M, n)
  alpn <- NULL
  dimnames(fred$ans) <- list(alpn, alpn, dimnames(m)[[1]])
  fred$a
}  # m2a






 a2m <- function(a, trim = FALSE) {

     
  if (is.matrix(a) && ncol(a) == nrow(a))
    a <- array(a, c(nrow(a), ncol(a), 1))
  if (!is.array(a))
    dim(a) <- c(1, 1, length(a))

  M <- dim(a)[1]
  if (M != dim(a)[2])
    stop("argument 'a' does not contain square matrices")
  n <- dim(a)[3]
  dimm.value <- dimm(M)
  index <- iam(NA, NA, M, both = TRUE, diag = TRUE)

  mat <- matrix(0, n, dimm.value)
  for (jay in seq(dimm.value))
    mat[, jay] <- a[index$row[jay], index$col[jay], ]

  if (trim)
    for (jay in dimm.value:1) {
      if (all(mat[, jay] == 0)) mat <- mat[, -jay] else break
    }

  mat
}  # a2m






 vindex <- function(M, row.arg = FALSE, col.arg = FALSE,
                    length.arg = M * (M + 1) / 2) {



  if ((row.arg + col.arg) != 1)
    stop("only one of row and col must be TRUE")
  if (M == 1) {
    ans <- 1
  } else {
    if (row.arg) {
      i1 <- matrix(1:M, M, M)
      ans <- c(i1[row(i1) + col(i1) <= (M + 1)])
    } else {
      i1 <- matrix(1:M, M, M)
      ans <- c(i1[row(i1) >= col(i1)])
    }
  }
  if (length.arg > length(ans))
    stop("argument 'length.arg' too big")
  rep_len(ans, length.arg)
}  # vindex






 wweights <-
  function(object, matrix.arg = TRUE, deriv.arg = FALSE,
           ignore.slot = FALSE, checkwz = TRUE) {






  if (length(wz <- object@weights) && !ignore.slot && !deriv.arg) {
    return(wz)
  }

  M <- object@misc$M  # Done below
  n <- object@misc$n  # Done below

  if (any(slotNames(object) == "extra")) {
    extra <- object@extra
    if (length(extra) == 1 && !length(names(extra))) {
      extra <- extra[[1]]
    }
  }
  mu <- object@fitted.values
  if (any(slotNames(object) == "predictors"))
    eta <- object@predictors
  mt <- terms(object)  # object@terms$terms; 20030811
  Hlist <- object@constraints
  new.coeffs <- object@coefficients
  if (any(slotNames(object) == "iter"))
    iter <- object@iter

  w <- rep_len(1, n)
  if (any(slotNames(object) == "prior.weights"))
    w <- object@prior.weights
  if (!length(w))
    w <- rep_len(1, n)

  x <- object@x
  if (!length(x))
    x <- model.matrixvlm(object, type = "lm")

  y <- object@y
  if (!length(y))
    y <- depvar(object)

    offset <- object@offset  # Could be cbind(0)
    if (all(dim(offset) == 1)) offset <- c(offset)  # 0 (not a matrix)


  X.vlm.save <- model.matrixvlm(object, type = "vlm")


  if (length(object@misc$form2)) {
    Xm2 <- object@Xm2
    if (!length(Xm2))
      Xm2 <- model.matrix(object, type = "lm2")
    Ym2 <- object@Ym2
  }




  if (any(slotNames(object) == "family")) {
    infos.list <- object@family@infos()
    if (length(infos.list))
      for (ii in names(infos.list)) {
        assign(ii, infos.list[[ii]])
      }
  }



  if (any(slotNames(object) == "control"))
    for (ii in names(object@control)) {
      assign(ii, object@control[[ii]])
    }

  if (length(object@misc))
    for (ii in names(object@misc)) {
      assign(ii, object@misc[[ii]])
    }




  if (any(slotNames(object) == "family")) {
    expr <- object@family@deriv
    deriv.mu <- eval(expr)
    if (!length(wz)) {
      expr <- object@family@weight
      wz <- eval(expr)


      if (M > 1)
        dimnames(wz) <- list(dimnames(wz)[[1]], NULL)  # Remove colnames
      wz <- if (matrix.arg) as.matrix(wz) else c(wz)
    }
    if (deriv.arg) list(deriv = deriv.mu, weights = wz) else wz
  } else {
    NULL
  }
}  # wweights






 pweights <- function(object, ...) {
  ans <- object@prior.weights
  if (length(ans)) {
    ans
  } else {
    temp <- object@y
    ans <- rep_len(1, nrow(temp))  # Assumed all equal and unity.
    names(ans) <- dimnames(temp)[[1]]
    ans
  }
}  # pweights









procVec <- function(vec, yn, Default) {








  if (anyNA(vec))
    stop("vec cannot contain any NAs")
  L <- length(vec)
  nvec <- names(vec)  # vec[""] undefined
  named <- length(nvec)  # FALSE for c(1,3)
  if (named) {
    index <- (1:L)[nvec == ""]
    default <- if (length(index)) vec[index] else Default
  } else {
    default <- vec
  }

  answer <- rep_len(default, length(yn))
  names(answer) <- yn
  if (named) {
    nvec2 <- nvec[nvec != ""]
    if (length(nvec2)) {
      if (any(!is.element(nvec2, yn)))
          stop("some names given which are superfluous")
      answer <- rep_len(NA_real_, length(yn))
      names(answer) <- yn
      answer[nvec2] <- vec[nvec2]
      answer[is.na(answer)] <-
        rep_len(default, sum(is.na(answer)))
    }
  }

  answer
}  # procVec



if (FALSE) {


}





 weightsvglm <-
  function(object, type = c("prior", "working"),
           matrix.arg = TRUE, ignore.slot = FALSE,
           deriv.arg = FALSE, ...) {
  weightsvlm(object, type = type, matrix.arg = matrix.arg,
              ignore.slot = ignore.slot,
              deriv.arg = deriv.arg, ...)
}



 weightsvlm <-
  function(object, type = c("prior", "working"),
           matrix.arg = TRUE, ignore.slot = FALSE,
           deriv.arg = FALSE, ...) {
  if (mode(type) != "character" && mode(type) != "name")
    type <- as.character(substitute(type))
  type <- match.arg(type, c("prior", "working"))[1]

  if (type == "working") {
    wweights(object = object,
             matrix.arg = matrix.arg, deriv.arg = deriv.arg,
             ignore.slot = ignore.slot, ...)
  } else {
    if (deriv.arg)
      stop("cannot set 'deriv = TRUE' when 'type=\"prior\"'")
    ans <- pweights(object)
    if (matrix.arg) as.matrix(ans) else c(ans)
  }
}  # weightsvlm



if (!isGeneric("weights"))
    setGeneric("weights", function(object, ...)
  standardGeneric("weights"))


setMethod("weights", "vlm",
         function(object, ...)
         weightsvlm(object, ...))


setMethod("weights", "vglm",
         function(object, ...)
         weightsvglm(object, ...))












qnupdate <- function(w, wzold, dderiv, deta, M, keeppd = TRUE,
                     trace = FALSE, reset = FALSE,
                     effpos=.Machine$double.eps^0.75) {


  if (M == 1) {
    dderiv <- cbind(dderiv)
    deta <- cbind(deta)
  }
  Bs <- mux22(t(wzold), deta, M = M,
              upper = FALSE, as.matrix = TRUE)  # n x M
  sBs <- c( (deta * Bs) %*% rep_len(1, M) )  # should have positive vals
  sy <- c( (dderiv * deta) %*% rep_len(1, M) )
  wznew <- wzold
  index <- iam(NA, NA, M = M, both = TRUE)
  index$row.index <- rep_len(index$row.index, ncol(wzold))
  index$col.index <- rep_len(index$col.index, ncol(wzold))
  updateThese <- if (keeppd) (sy > effpos) else rep_len(TRUE, length(sy))
  if (!keeppd || any(updateThese)) {
    wznew[updateThese,] <- wznew[updateThese,] -
        Bs[updateThese,index$row] *
        Bs[updateThese,index$col] / sBs[updateThese] +
        dderiv[updateThese,index$row] *
        dderiv[updateThese,index$col] / sy[updateThese]
    notupdated <- sum(!updateThese)
    if (notupdated && trace)
      cat(notupdated,
          "weight matrices not updated out of", length(sy), "\n")
  } else {
    warning("no BFGS quasi-Newton update made at all")
    cat("no BFGS quasi-Newton update made at all\n")
    flush.console()
  }
  wznew
}  # qnupdate










mbesselI0 <-
  function(x, deriv.arg = 0) {
  if (!is.Numeric(deriv.arg, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE) &&
      deriv.arg != 0)
    stop("argument 'deriv.arg' must be a single non-negative integer")
  if (!(deriv.arg == 0 || deriv.arg == 1 || deriv.arg == 2))
    stop("argument 'deriv' must be 0, 1, or 2")
  if (!is.Numeric(x))
    stop("bad input for argument 'x'")
  nn <- length(x)

  if (FALSE) {
    }

    ans <- matrix(NA_real_, nrow = nn, ncol = deriv.arg+1)
    ans[, 1] <- besselI(x, nu = 0)
    if (deriv.arg>=1) ans[, 2] <- besselI(x, nu = 1)
    if (deriv.arg>=2) ans[, 3] <- ans[,1] - ans[,2] / x
    ans
}  # mbesselI0






VGAM.matrix.norm <- function(A, power = 2, suppressWarning = FALSE) {
  if ((nrow(A) != ncol(A)) && !suppressWarning)
    warning("norms should be calculated for square matrices; ",
            "'A' is not square")
  if (power == "F") {
    sqrt(sum(A^2))
  } else if (power == 1) {
    max(colSums(abs(A)))
  } else if (power == 2) {
    sqrt(max(eigen(t(A) %*% A, symmetric = TRUE)$value))
  } else if (!is.finite(power)) {
    max(colSums(abs(A)))
  } else {
    stop("argument 'power' not recognized")
  }
}  # VGAM.matrix.norm






rmfromVGAMenv <- function(varnames, prefix = "") {
  evarnames <- paste(prefix, varnames, sep = "")
  for (ii in evarnames) {
    mytext1 <- "exists(x = ii, envir = VGAMenv)"
    myexp1 <- parse(text = mytext1)
    is.there <- eval(myexp1)
    if (is.there) {
      rm(list = ii, envir = VGAMenv)
    }
  }
}  # rmfromVGAMenv






existsinVGAMenv <- function(varnames, prefix = "") {
  evarnames <- paste(prefix, varnames, sep = "")
  ans <- NULL
  for (ii in evarnames) {
    mytext1 <- "exists(x = ii, envir = VGAMenv)"
    myexp1 <- parse(text = mytext1)
    is.there <- eval(myexp1)
    ans <- c(ans, is.there)
  }
  ans
}  # existsinVGAMenv






assign2VGAMenv <- function(varnames, mylist, prefix = "") {
  evarnames <- paste(prefix, varnames, sep = "")
  for (ii in seq_along(varnames)) {
    assign(evarnames[ii], mylist[[(varnames[ii])]],
           envir = VGAMenv)
  }
}  # assign2VGAMenv






getfromVGAMenv <- function(varname, prefix = "") {
  varname <- paste(prefix, varname, sep = "")
  if (length(varname) > 1)
    stop("'varname' must be of length 1")
  get(varname, envir = VGAMenv)
}






lerch <- function(x, s, v, tolerance = 1.0e-10, iter = 100) {
  if (!is.Numeric(x) || !is.Numeric(s) || !is.Numeric(v))
    stop("bad input in 'x', 's', and/or 'v'")
  if (is.complex(c(x,s,v)))
    stop("complex arguments not allowed in 'x', 's' and 'v'")
  if (!is.Numeric(tolerance, length.arg = 1, positive = TRUE) ||
      tolerance > 0.01)
    stop("bad input for argument 'tolerance'")
  if (!is.Numeric(iter, length.arg = 1,
                  integer.valued = TRUE, positive = TRUE))
    stop("bad input for argument 'iter'")

  L <- max(length(x), length(s), length(v))
  x <- rep_len(x, L)
  s <- rep_len(s, L)
  v <- rep_len(v, L)
  xok <- abs(x) < 1 & !(v <= 0 & v == round(v))
  x[!xok] <- 0  # Fix this later

  ans <- .C("lerchphi123",
           err = integer(L), as.integer(L),
           as.double(x), as.double(s), as.double(v),
           acc=as.double(tolerance), result=double(L),
           as.integer(iter))

  ifelse(ans$err == 0 & xok , ans$result, NA)
}  # lerch




negzero.expression.VGAM <- expression({










  if (length(dotzero) == 1 && (dotzero == "" || is.na(dotzero)))
    dotzero <- NULL



  if (is.character(dotzero)) {




  which.numeric.all <- NULL
  for (ii in seq_along(dotzero)) {
    which.ones <-
      grep(dotzero[ii], predictors.names, fixed = TRUE)
    if (length(which.ones)) {
      which.numeric.all <- c(which.numeric.all, which.ones)
    } else {
      warning("some values of argument 'zero' are unmatched. ",
              "Ignoring them")
    }
  }
  which.numeric <- unique(sort(which.numeric.all))

  if (!length(which.numeric)) {
    warning("No values of argument 'zero' were matched.")
    which.numeric <- NULL
  } else if (length(which.numeric.all) > length(which.numeric)) {
    warning("There were redundant values of argument 'zero'.")
  }

    dotzero <- which.numeric
  }



  posdotzero <-  dotzero[dotzero > 0]
  negdotzero <-  dotzero[dotzero < 0]

  zneg.index <- if (length(negdotzero)) {

    if (!is.Numeric(-negdotzero, positive = TRUE,
                    integer.valued = TRUE) ||
        max(-negdotzero) > M1)
        stop("bad input for argument 'zero'")

    bigUniqInt <- 1080
    zneg.index <- rep(0:bigUniqInt, rep(length(negdotzero),
                      1 + bigUniqInt)) * M1 + abs(negdotzero)
    sort(intersect(zneg.index, 1:M))
  } else {
    NULL
  }

  zpos.index <- if (length(posdotzero)) posdotzero else NULL
  z.Index <- if (!length(dotzero))
               NULL else
               unique(sort(c(zneg.index, zpos.index)))

  constraints <- cm.zero.VGAM(constraints, x = x, z.Index, M = M,
                              predictors.names = predictors.names,
                              M1 = M1)
})  # negzero.expression.VGAM






is.empty.list <- function(mylist) {
  is.list(mylist) &&
  length(unlist(mylist)) == 0
}









  interleave.VGAM  <- function(.M, M1, inverse = FALSE) {
  if (inverse) {
    NRs <- (.M)/M1
    if (round(NRs) != NRs)
      stop("Incompatible number of parameters")
    c(matrix(1:(.M), nrow = NRs, byrow = TRUE))
  } else {
    c(matrix(1:(.M), nrow = M1, byrow = TRUE))
  }
}  # interleave.VGAM








interleave.cmat <- function(cmat1, cmat2) {
  ncol1 <- ncol(cmat1)
  ncol2 <- ncol(cmat2)
  if (ncol1 == 1) {
    return(cbind(cmat1, cmat2))
  } else {  # ncol1 > 1
    if (ncol2 == 1) {
      return(cbind(cmat1[, 1], cmat2, cmat1[, -1]))
    } else
    if (ncol1 != ncol2) {
      warning("this function is confused. ",
              "Returning cbind(cmat1, cmat2)")
      return(cbind(cmat1[, 1], cmat2, cmat1[, -1]))
    } else {  # ncol1 == ncol2 and both are > 1.
      kronecker(cmat1, cbind(1, 0)) +
      kronecker(cmat2, cbind(0, 1))
    }
  }
}  # interleave.cmat






w.wz.merge <- function(w, wz, n, M, ndepy,
                       intercept.only = FALSE) {





  wz <- as.matrix(wz)

  if (ndepy == 1)
    return( c(w) * wz)


  if (intercept.only)
    warning("yettodo: support intercept.only == TRUE")

  if (NCOL(w) > ndepy)
    stop("number of columns of 'w' exceeds number of responses")

  w  <- matrix(w, n, ndepy)
  w.rep <- matrix(0, n, ncol(wz))
  M1 <- M / ndepy
  all.indices <- iam(NA, NA, M = M, both = TRUE)



  if (FALSE)
  for (ii in 1:ncol(wz)) {

    if ((ind1 <- ceiling(all.indices$row[ii] / M1)) ==
                 ceiling(all.indices$col[ii] / M1)) {
      w.rep[, ii] <- w[, ind1]
    }


  }  # ii


  res.Ind1 <- ceiling(all.indices$row.index / M1)
  Ind1 <- res.Ind1 == ceiling(all.indices$col.index / M1)

  LLLL <- min(ncol(wz), length(Ind1))
  Ind1 <- Ind1[1:LLLL]
  res.Ind1 <- res.Ind1[1:LLLL]

  for (ii in 1:ndepy) {
    sub.ind1 <- (1:LLLL)[Ind1 & (res.Ind1 == ii)]
    w.rep[, sub.ind1] <- w[, ii]
  }  # ii

  w.rep * wz
}  # w.wz.merge





w.y.check <- function(w, y,
                      ncol.w.max = 1, ncol.y.max = 1,
                      ncol.w.min = 1, ncol.y.min = 1,
                      out.wy = FALSE,
                      colsyperw = 1,
                      maximize = FALSE,
                      Is.integer.y = FALSE,
                      Is.positive.y = FALSE,
                      Is.nonnegative.y = FALSE,
                      prefix.w = "PriorWeight",
                      prefix.y = "Response") {



  if (!is.matrix(w))
    w <- as.matrix(w)
  if (!is.matrix(y))
    y <- as.matrix(y)
  n.lm <- nrow(y)
  rn.w <- rownames(w)
  rn.y <- rownames(y)
  cn.w <- colnames(w)
  cn.y <- colnames(y)


  if (Is.integer.y && any(y != round(y)))
    stop("response variable 'y' must be integer-valued")
  if (Is.positive.y && any(y <= 0))
    stop("response variable 'y' must be positive-valued")
  if (Is.nonnegative.y && any(y < 0))
    stop("response variable 'y' must be 0 or positive-valued")

  if (nrow(w) != n.lm)
    stop("nrow(w) should be equal to nrow(y)")

  if (ncol(w) > ncol.w.max)
    stop("prior-weight variable 'w' has too many columns")
  if (ncol(y) > ncol.y.max)
    stop("response variable 'y' has too many columns; ",
         "only ", ncol.y.max, " allowed")

  if (ncol(w) < ncol.w.min)
    stop("prior-weight variable 'w' has too few columns")
  if (ncol(y) < ncol.y.min)
    stop("response variable 'y' has too few columns; ",
         "at least ", ncol.y.max, " needed")

  if (min(w) <= 0)
    stop("prior-weight variable 'w' must contain positive ",
         "values only")

  if (is.numeric(colsyperw) && ncol(y) %% colsyperw != 0)
    stop("number of columns of the response variable 'y' is not ",
         "a multiple of ", colsyperw)


  if (maximize) {
    Ncol.max.w <- max(ncol(w), ncol(y) / colsyperw)
    Ncol.max.y <- max(ncol(y), ncol(w) * colsyperw)
  } else {
    Ncol.max.w <- ncol(w)
    Ncol.max.y <- ncol(y)
  }

  if (out.wy && ncol(w) < Ncol.max.w) {
    nblanks <- sum(cn.w == "")
    if (nblanks > 0)
      cn.w[cn.w == ""] <- param.names(prefix.w, nblanks)
    if (length(cn.w) < Ncol.max.w)
      cn.w <- c(cn.w, paste(prefix.w, (length(cn.w)+1):Ncol.max.w,
                            sep = ""))
    w <- matrix(w, n.lm, Ncol.max.w, dimnames = list(rn.w, cn.w))
  }
  if (out.wy && ncol(y) < Ncol.max.y) {
    nblanks <- sum(cn.y == "")
    if (nblanks > 0)
      cn.y[cn.y == ""] <- param.names(prefix.y, nblanks)
    if (length(cn.y) < Ncol.max.y)
      cn.y <- c(cn.y, paste(prefix.y, (length(cn.y)+1):Ncol.max.y,
                            sep = ""))
    y <- matrix(y, n.lm, Ncol.max.y, dimnames = list(rn.y, cn.y))
  }

  list(w = if (out.wy) w else NULL,
       y = if (out.wy) y else NULL)
}  # w.y.check






arwz2wz <-
    function(arwz, M = 1, M1 = 1, rm.trailing.cols = TRUE,
             full.arg = FALSE) {



  if (length(dim.arwz <- dim(arwz)) != 3)
    stop("dimension of 'arwz' should be of length 3")
  n       <- dim.arwz[1]
  ndepy   <- dim.arwz[2]
  dim.val <- dim.arwz[3]

  if (ndepy == 1) {
    dim(arwz) <- c(n, dim.val)
    return(arwz)
  }

  wz <- matrix(0.0, n,
               if (full.arg) M*(M+1)/2 else sum(M:(M-M1+1)))
  ind1 <- iam(NA, NA, M = M1, both = TRUE, diag = TRUE)
  len.ind1 <- dim.val # length(ind1$col.index)

  for (ii in 1:ndepy) {
    for (jlocal in 1:len.ind1) {
      wz[, iam(M1 * (ii - 1) + ind1$row[jlocal],
               M1 * (ii - 1) + ind1$col[jlocal],
               M = M)] <- arwz[, ii, jlocal]
    }
  }

  if (rm.trailing.cols && !full.arg) {
    colind <- ncol(wz)
    while (all(wz[, colind] == 0))
      colind <- colind - 1
    if (colind < ncol(wz))
      wz <- wz[, 1:colind, drop = FALSE]
  }

  wz
}  # arwz2wz






 wz.merge <-
    function(wz1, wz2, M1, M2, rm.trailing.cols = TRUE) {


  if (!is.matrix(wz1))
    wz1 <- cbind(wz1)
  if (!is.matrix(wz2))
    wz2 <- cbind(wz2)
  M <- M1 + M2
  wz <- matrix(0.0, nrow(wz1), M*(M+1)/2)

  for (ilocal in 1:M1)
    for (jlocal in ilocal:M1)
      if (iam(ilocal, jlocal, M = M1) <= ncol(wz1)) {
         wz[, iam(ilocal, jlocal, M = M)] <-
        wz1[, iam(ilocal, jlocal, M = M1)]
      }

  for (ilocal in 1:M2)
    for (jlocal in ilocal:M2)
      if (iam(ilocal, jlocal, M = M2) <= ncol(wz2)) {
         wz[, iam(M1+ilocal, M1+jlocal, M = M)] <-
        wz2[, iam(   ilocal,    jlocal, M = M2)]
      }

  if (rm.trailing.cols) {
    colind <- ncol(wz)
    while (all(wz[, colind] == 0))
      colind <- colind - 1
    if (colind < ncol(wz))
      wz <- wz[, 1:colind, drop = FALSE]
  }

  wz
}  # wz.merge







param.names <-
    function(string, S = 1, skip1 = FALSE, sep = "") {
  if (skip1) {
    if (S == 1) string else paste(string, 1:S, sep = sep)
  } else {
    paste(string, 1:S, sep = sep)
  }
}  # param.names







vweighted.mean.default <-
    function (x, w, ..., na.rm = FALSE) {
  tmp5 <- w.y.check(w = w, y = x, ncol.w.max = Inf,
                    ncol.y.max = Inf,
                    out.wy = TRUE,
                    colsyperw = 1,
                    maximize = TRUE,
                    Is.integer.y = FALSE,
                    Is.positive.y = FALSE,
                    Is.nonnegative.y = FALSE,
                    prefix.w = "PriorWeight",
                    prefix.y = "Response")

  x <- tmp5$y
  w <- tmp5$w

  ans <- numeric(ncol(w))
  for (ii in 1:ncol(w))
      ans[ii] <- weighted.mean(x[, ii], w = w[, ii], ...,
                               na.rm = na.rm)
  ans
}  # vweighted.mean.default






familyname.vlm <- function(object, all = FALSE, ...) {
  ans <- object@family@vfamily
  if (all) ans else ans[1]
}


familyname.vglmff <- function(object, all = FALSE, ...) {
  ans <- object@vfamily
  if (all) ans else ans[1]
}



if (!isGeneric("familyname"))
    setGeneric("familyname",
  function(object, ...) standardGeneric("familyname"))


setMethod("familyname", "vglmff",
         function(object, ...)
         familyname.vglmff(object, ...))



setMethod("familyname", "vlm",
         function(object, ...)
         familyname.vlm(object, ...))










bisection.basic <-
  function(f, a, b, tol = 1e-9, nmax = NULL, ...) {




      

  if (any(is.infinite(b))) {
    warning("replacing 'b' values of Inf by a large value")
    b[is.infinite(b)] <- .Machine$double.xmax / 4
  }

 


  if (is.null(nmax)) {
    nmax <- round(log2(max(b - a)) - log2(min(tol))) + 4
    if (!is.finite(nmax))
      nmax <- log2(.Machine$double.xmax) - 5
  }
  signtest <- (sign(f(a, ...)) * sign(f(b, ...)) <= 0)

  allsign <- all(signtest, na.rm = TRUE)
  if (!allsign || any(is.na(signtest))){
    warning("roots do not exist between 'a' and 'b'. ",
            "Some answers may be misleading.")
  }
      


  N <- 1
  while (N <= nmax) {
    mid <- (a + b) / 2
    save.f <- f(mid, ...)
    if (all(save.f == 0 | (b - a)/2 < tol)) {
      return(mid)
    }
    N <- N + 1
    vecTF <- sign(save.f) == sign(f(a, ...))
    a[ vecTF] <- mid[ vecTF]
    b[!vecTF] <- mid[!vecTF]
  }

  warning("did not coverge. Returning final root")
  mid
}  # bisection.basic





retain.col <- function(mat, coln
                      ) {
  if (is.matrix(mat))  # && exclude
    mat[, -coln] <- 0
  mat
}  # retain.col






which.etas <-
  function(object, kay = 1) {

  cmat <- constraints(object, matrix = TRUE)
  if (ncol(cmat) < kay)
    stop("value of argument 'kay' is too large")
  which(cmat[, kay] != 0)
}  # which.etas





which.xij <-
  function(object, ...) {

  cmat <- constraints(object, matrix = TRUE)
  ans <- rep(FALSE, NCOL(cmat))
  names(ans) <- colnames(cmat)
  xij <- object@control$xij
  if (!is.null(xij)) {
    for (ii in 1:length(xij)) {
      responsevar <- all.vars(xij[[ii]])[1]
      ans[responsevar] <- TRUE  # Assumes NCOL(cmat) == 1
    }
  }
  ans
}  # which.xij







attr.assign.x.vglm <- function(object) {
  x.lm.vglm <- NULL
  x.lm.vglm <- try(model.matrix(formula(object),
                                data = get(object@misc$dataname)), TRUE)
  if (!length(x.lm.vglm))
    x.lm.vglm <- try(model.matrix(formula(object)), TRUE)
  if (!length(x.lm.vglm))
    stop("cannot create the 'lm'-type model matrix from formula(object)")


  attr.x.lm.vglm <- attr(x.lm.vglm, "assign")

  clist <- constraints(object, type = "term")  # type = "lm" is default
  ncols.cm <- unlist(lapply(clist, function(cm) ncol(cm)))

  icounter <- 1
  attr.x.vglm <- NULL
  for (ii in min(attr.x.lm.vglm):max(attr.x.lm.vglm)) {
    attr.x.vglm <- c(attr.x.vglm,
                     rep(ii, sum(attr.x.lm.vglm == ii) *
                             ncols.cm[icounter]))
    icounter <- icounter + 1
  }

  attr.x.vglm
}  # attr.assign.x.vglm







muxtXAX <- function(A, X, M) {

  if ((n <- nrow(X)) != nrow(A))
    stop("arguments 'X' and 'A' are not conformable")

  ans1 <- array(0, c(n, M, M))
  for (eye in 1:M)
    for (jay in 1:M)
      for (kay in 1:M)
        ans1[, eye, jay] <- ans1[, eye, jay] +
          A[, iam(eye, kay, M)] * X[, iam(kay, jay, M)]

  ans2 <- matrix(0, n, dimm(M))
  for (eye in 1:M)
    for (jay in eye:M)
      for (kay in 1:M)
        ans2[, iam(eye, jay, M)] <- ans2[, iam(eye, jay, M)] +
        X[, iam(kay, eye, M)] * ans1[, kay, jay]
  ans2
}  # muxtXAX











 trim.constraints <-
   function(object, sig.level = 0.05,
            max.num = Inf,
            intercepts = TRUE,   # intercepts constraint matrices
            ...) {





  if (!is(object, "vlm"))
    stop("argument 'object' should inherit from the 'vglm' class")


  if (!is.finite(sig.level) || !is.numeric(sig.level) ||
      min(sig.level) < 0 || max(sig.level) > 1)
    stop("bad input for argument 'sig.level'")

  csfit <- coef(summary(object, HDEtest = FALSE))
  X.small <- model.matrix(object, type = "lm")
  X.assign <- attr(X.small, "assign")
  cmlist0 <-
  cmlist2 <-
  cmlist1 <- constraints(object, type = "term")
  ncolHlist <- sapply(cmlist1, ncol)
  colsperterm <- sapply(X.assign, length)
  endpts <- cumsum(colsperterm * ncolHlist)
  sig.level <- rep_len(sig.level, max(endpts))


  pvs.vec <- csfit[, "Pr(>|z|)"]
  if (any(!is.finite(pvs.vec)))
    stop("cannot handle any NAs, etc. as p-values")
  if (is.finite(max.num)) {
    if (!is.Numeric(max.num, positive = TRUE, integer.valued = TRUE))
      stop("bad input for argument 'max.num'")
    spvs.vec <- sort(pvs.vec, decreasing = TRUE)
    one.more <- as.numeric(max.num < length(spvs.vec))
    use.sig.level <- mean(spvs.vec[max.num + (0:one.more)])
    sig.level[sig.level < use.sig.level] <- use.sig.level
  }

  for (kay in rev(seq(length(cmlist1)))) {
    cm.kay0 <- cm.kay <- cmlist1[[kay]]
    for (cay in rev(seq(ncol(cm.kay0)))) {
      cptr <- (if (kay == 1) 0 else endpts[kay - 1]) +
              (ncolHlist[kay]) * seq(colsperterm[kay]) +
              cay - ncol(cm.kay0)
      if (all(pvs.vec[cptr] > sig.level[cptr]))
        cm.kay <- cm.kay[, -cay, drop = FALSE]
    }  # cay
    cmlist2[[kay]] <- if (length(cm.kay)) cm.kay else NULL
  }




  if (!intercepts && length(cmlist0[["(Intercept)"]])) {
    cmlist2[["(Intercept)"]] <- cmlist0[["(Intercept)"]]
  }



  cmlist2
}  # trim.constraints