.packageName <- "geepack"
geese <- function(formula = formula(data),
                  sformula = ~ 1,
                  id, waves = NULL,
                  data = parent.frame(), subset = NULL, na.action = na.omit,
                  contrasts = NULL, weights = NULL,
                  ## zcor is design matrix for alpha,
                  ## corp is known paratemers to correlation coef. rho
                  zcor = NULL, corp = NULL,
                  ## zsca is constructed from sformula
                  ## control parameters
                  control = geese.control(...),
                  ## param 
                  b = NULL, alpha = NULL, gm = NULL,
                  ## geestr
                  family = gaussian(),
                  mean.link = NULL,
                  variance = NULL,
                  cor.link = "identity",
                  sca.link = "identity",
                  link.same = TRUE,
                  scale.fix = FALSE, scale.value = 1.0,
                  ## corr
                  corstr = "independence",
                  ...) {
  scall <- match.call()
  mnames <- c("", "formula", "data", "offset", "weights", "subset", "na.action", "id", "waves", "corp")
  cnames <- names(scall)
  cnames <- cnames[match(mnames,cnames,0)]
  mcall <- scall[cnames]
  if (is.null(mcall$id)) mcall$id <- as.name("id")
  mcall[[1]] <- as.name("model.frame")
  m <- eval(mcall, parent.frame())

  y <- model.extract(m, response)
  if (is.null(dim(y))) N <- length(y) else N <- dim(y)[1]
  mterms <- attr(m, "terms")
  x <- model.matrix(mterms, m, contrasts)
  offset <- model.extract(m, offset)
  if (is.null(offset)) offset <- rep(0, N)
  w <- model.extract(m, weights)
  if (is.null(w)) w <- rep(1, N)
  id <- model.extract(m, id)
  waves <- model.extract(m, waves)
  corp <- model.extract(m, corp)
  if (is.null(id)) stop("id variable not found.")

  ## setting up the scale model;
  ## borrowed idea from S+ function dglm by Gordon Smyth
  mcall$formula <- formula
  mcall$formula[3] <- switch(match(length(sformula), c(0,2,3)),
                             1, sformula[2], sformula[3])
  m <- eval(mcall, parent.frame())
  terms <- attr(m, "terms")
  zsca <- model.matrix(terms, m, contrasts)
  soffset <- model.extract(m, offset)
  if (is.null(soffset)) soffset <- rep(0, N)  
 
  if (is.character(family)) family <- get(family)
  if (is.function(family))  family <- family()
  ans <- geese.fit(x, y, id, offset, soffset, w,
                   waves, zsca, zcor, corp, 
                   control,
                   b, alpha, gm,
                   family, mean.link, variance, cor.link, sca.link,
                   link.same, scale.fix, scale.value, 
                   corstr, ...)
  ans <- c(ans, list(call=scall, formula=formula)) 
  class(ans) <- "geese"
  ans
}

geese.fit <- function(x, y, id,
                      offset=rep(0,N), soffset=rep(0,N), weights=rep(1,N),
                      waves = NULL, zsca = matrix(1,N,1),
                      zcor = NULL, corp = NULL,
                      control = geese.control(...),
                      ## param 
                      b = NULL, alpha = NULL, gm = NULL,
                      ## geestr
                      family = gaussian(),
                      mean.link = NULL,
                      variance = NULL,
                      cor.link = "identity",
                      sca.link = "identity",
                      link.same = TRUE,
                      scale.fix = FALSE, scale.value = 1.0,
                      ## corr
                      corstr = "independence", ...) {
  N <- length(id)
  ##clusz <- unlist(lapply(split(id, id), length))
  clusnew <- c(which(diff(as.numeric(id)) != 0), length(id))
  clusz <- c(clusnew[1], diff(clusnew))
  maxclsz <- max(clusz)
  if (is.null(waves)) waves <- unlist(sapply(clusz, function(x) 1:x))
  waves <- as.integer(waves)

  LINKS <- c("identity", "logit", "probit", "cloglog", "log", "inverse", "fisherz", "lwybc2", "lwylog")
  VARIANCES <- c("gaussian", "binomial", "poisson", "Gamma") ## quasi is not supported yet

  if (is.null(mean.link)) mean.link <- family$link
  if (is.null(variance)) variance <- family$family
  mean.link.v <- pmatch(mean.link, LINKS, -1, TRUE)
  cor.link.v <- pmatch(cor.link, LINKS, -1, TRUE)
  sca.link.v <- pmatch(sca.link, LINKS, -1, TRUE)
  variance.v <- pmatch(variance, VARIANCES, -1, TRUE)
  if (any(mean.link.v == -1)) stop("mean.link invalid.")
  if (any(cor.link.v == -1)) stop("cor.link invalid.")
  if (any(sca.link.v == -1)) stop("sca.link invalid.")
  if (any(variance.v == -1)) stop("variance invalid.")
  if (length(mean.link.v) != length(variance.v))
    stop("mean.link and variance not same length.")
  if (length(mean.link.v) != length(sca.link.v))
    stop("mean.link and sca.link not same lehgnt.")
      
  if (length(id) != length(y)) stop("id and y not same length.")
  if (length(offset) != length(y)) stop("offset and y not same length")
  if (length(soffset) != length(y)) stop("sca.offset and y not same length")
  if (nrow(zsca) != length(y)) stop("nrow(zsca) and length(y) not match")
  
  if (link.same) linkwaves <- rep(1, N)
  else {
    if (max(waves) != maxclsz) stop("maximum waves and maximum cluster size not equal")
    if (length(mean.link.v) != maxclsz) stop("length of mean.link not equal to the maximum cluster size.")
    linkwaves <- waves
  }
  linkwaves <- as.integer(linkwaves)
  geestr <- list(length(mean.link.v), as.integer(mean.link.v),
                 as.integer(variance.v), as.integer(sca.link.v),
                 as.integer(cor.link.v), as.integer(scale.fix))

  CORSTRS <- c("independence", "exchangeable", "ar1", "unstructured", "userdefined")
  corstrv <- pmatch(corstr, CORSTRS, -1)
  if (corstrv == -1) stop("invalid corstr.")
  corr <- list(as.integer(corstrv), maxclsz)
  
  if (is.null(zcor)) {
    if (corstrv == 5) stop("need zcor matrix for userdefined corstr.") 
    else zcor <- genZcor(clusz, waves, corstrv)
  }
  else {
    if (!is.matrix(zcor)) zcor <- as.matrix(zcor)
    if (corstrv >= 4 && nrow(zcor) != sum(clusz * (clusz - 1) / 2)) stop("nrow(zcor) need to be equal sum(clusz * (clusz - 1) / 2) for unstructured or userdefined corstr.")
    if (corstrv %in% c(2,3) && nrow(zcor) != length(clusz)) stop("nrow(zcor) need to be equal to the number of clusters for exchangeable or ar1 corstr.")
  }
  if (!is.matrix(zcor)) zcor <- as.matrix(zcor)
  if (is.null(corp)) corp <- as.double(waves)

  p <- ncol(x)
  q <- ncol(zcor)
  r <- ncol(zsca)
  
  ## Initial values setup
  fit0 <- glm.fit(x, y, weights=weights, offset=offset, family=family)
  if (is.null(b)){
    ##b <- rep(1,p)
    b <- fit0$coef
  }
  if (is.null(alpha)) alpha <- rep(0,q)
  if (is.null(gm)) {
    ##gm <- rep(scale.value, r)
    qlf <- quasi(LINKS[sca.link.v])$linkfun
    pr2 <- (residuals.glm(fit0, type="pearson")) ^ 2
    gm <- lm.fit(zsca, qlf(pr2), offset = soffset)$coef
  }
  param <- list(b, alpha, gm)

  ans <- .Call("gee_rap", y, x, offset, soffset, weights,
               linkwaves, zsca, zcor, corp,
               clusz, geestr, corr, param, control, PACKAGE = "geepack")
  names(ans) <- c("beta", "alpha", "gamma", "vbeta", "valpha", "vgamma",
                  "vbeta.naiv", "valpha.naiv", "valpha.stab",
                  "vbeta.ajs", "valpha.ajs", "vgamma.ajs",
                  "vbeta.j1s", "valpha.j1s", "vgamma.j1s",
                  "vbeta.fij", "valpha.fij", "vgamma.fij",
                  "error")
  ans$xnames <- dimnames(x)[[2]]
  ans$zsca.names <- dimnames(zsca)[[2]]
  ans$zcor.names <- dimnames(zcor)[[2]]
  if (is.null(ans$zcor.names)) ans$zcor.names = paste("alpha", 1:ncol(zcor), sep=":")
  names(ans$beta) <- ans$xnames
  names(ans$gamma) <- ans$zsca.names
  if (length(ans$alpha) > 0)  names(ans$alpha) <- ans$zcor.names

  ans <- c(ans, list(clusz=clusz, control=control,
                     model=list(mean.link=mean.link,
                       variance=variance, sca.link=sca.link,
                       cor.link=cor.link, corstr=corstr, scale.fix=scale.fix)))
  ans
}

geese.control <- function (epsilon = 1e-04, maxit = 25, trace = FALSE,
                           scale.fix = FALSE, jack = FALSE,
                           j1s = FALSE, fij = FALSE) {
  if (!is.numeric(epsilon) || epsilon <= 0) 
    stop("value of epsilon must be > 0")
  if (!is.numeric(maxit) || maxit <= 0) 
    stop("maximum number of iterations must be > 0")
  list(trace = as.integer(trace),
       jack = as.integer(jack), j1s = as.integer(j1s), fij = as.integer(fij),
       maxit = as.integer(maxit), epsilon = epsilon)
}
crossutri <- function(wave) {
  n <- length(wave)
  if (n == 1) return(NULL)
  ans <- rep(0, n*(n-1)/2)
  k <- 1
  for (i in 1:(n-1))
    for (j in (i+1):n) {
      ans[k] <- paste(wave[i], wave[j], sep=":")
      k <- k + 1
    }
  ans
}

genZcor <- function(clusz, waves, corstrv) {
  if (corstrv == 1) return (matrix(0,0,0))
  crs <- clusz * (clusz - 1) / 2
  if (corstrv == 2 || corstrv == 3) {
    ans <-  matrix(1, length(clusz), 1)
    ##ans <-  matrix(1, sum(crs), 1)
    colnames(ans) <- c("alpha")
  }
  else {
    id <- rep(1:length(clusz), clusz)
    z <- as.factor(unlist(lapply(split(waves, id), crossutri)))
    ans <- model.matrix(~z - 1)
    znames <- paste("alpha", unlist(crossutri(1:max(clusz))), sep=".")
    ##dimnames(ans) <- list(1:sum(crs), znames)
    colnames(ans) <- znames
  }
  ans
}


genZodds <- function(clusz, waves, corstrv, ncat) {
  if (corstrv == 1) return (matrix(0,0,0))
  crs <- clusz * (clusz - 1) / 2
  c2 <- ncat * ncat
  if (corstrv == 2 | corstrv == 3) {
    ans <- matrix(1, sum(crs) * c2, 1)
    colnames(ans) <- c("alpha")
  }
  else {
    id <- rep(1:length(clusz), clusz)
    z <- as.factor(unlist(lapply(split(waves, id), crossutri)))
    z <- model.matrix(~z - 1)
    ind <- gl(sum(crs), c2)
    ans <- z[ind,]
    colnames(ans) <- paste("alpha", 1:dim(ans)[2], sep=".")
  }
  ans
}
ordgee <- function(formula = formula(data), ooffset = NULL,
                   id, waves = NULL,
                   data=parent.frame, subset=NULL, na.action=na.omit,
                   contrasts=NULL, weights=NULL,
                   z=NULL, ##family=binomial(),
                   mean.link="logit",
                   corstr="independence",
                   control=geese.control(...),
                   b=NA, alpha=NA,
                   scale.fix=TRUE, scale.val=1,
                   int.const=TRUE, rev=FALSE, ##rev TRUE for coding in HZ 1996.
                   ...) {
### y is sum(n_i) * c x 1
### x is sum(n_i) * c x (p + c)
  scall <- match.call()
  mnames <- c("", "formula", "data", "offset", "weights", "subset", "id", "waves")
  cnames <- names(scall)
  cnames <- cnames[match(mnames,cnames,0)]
  mcall <- scall[cnames]
  if (is.null(mcall$id)) mcall$id <- as.name("id")
  mcall[[1]] <- as.name("model.frame")
  m <- eval(mcall, parent.frame())

  id <- model.extract(m, id)
##  N <- length(unique(id))
  clusz <- unlist(lapply(split(id, id), length))
  maxclsz <- max(clusz)
  if (is.null(waves)) waves <- unlist(sapply(clusz, function(x) 1:x))
  else waves <- model.extract(m, waves)
#   if (is.na(b)){
#    foo <- polr(formula, data, ...)
#    b <- c(foo$zeta, foo$coef)
#   }

  y <- model.extract(m, response)
  if (length(y) != length(id)) stop("response and id are not of the same length.")
  lev <- levels(y)
  nlev <- length(lev)
  ncat <- nlev - 1
  y <- unclass(y)
  Y <- rep(y, rep(ncat, sum(clusz)))
  if (rev) Y <- as.double(Y <= rep(1:ncat, sum(clusz)))
  else Y <- as.double(Y > rep(1:ncat, sum(clusz)))
  
  mterms <- attr(m, "terms")
  x <- model.matrix(mterms, m, contrasts)
  xvars <- as.character(attr(mterms, "variables"))[-1]
  if ((yvar <- attr(mterms, "response")) > 0) 
    xvars <- xvars[-yvar]
  xlev <- if (length(xvars) > 0) {
    xlev <- lapply(m[xvars], levels)
    xlev[!sapply(xlev, is.null)]
  }
  xint <- match("(Intercept)", colnames(x), nomatch = 0)
  n <- nrow(x)
  pc <- ncol(x)
  if (xint > 0) {
    x <- x[, -xint, drop = FALSE]
    pc <- pc - 1
  }
  else warning("an intercept is needed and assumed")
  ind <- gl(sum(clusz), ncat)
  x <- x[ind,, drop=FALSE]
    
  if (int.const) {
    xc <- matrix(diag(ncat), sum(clusz) * ncat, ncat, byrow=TRUE)
    colnames(xc) <- paste("Inter", lev[1:ncat], sep=":")
  }
  else {
    foo <- sapply(waves,
                  function(x, maxclsz, ncat) {
                    bar <- matrix(0, maxclsz*ncat, ncat)
                    bar[(x-1)*ncat + 1:ncat,] <- diag(ncat)
                    bar
                  }, maxclsz=maxclsz, ncat=ncat)
    xc <- matrix(unlist(foo), ncol=maxclsz*ncat, byrow=TRUE)
    colnames(xc) <- paste("Inter", paste(rep(1:maxclsz, rep(ncat, maxclsz)), rep(lev[1:ncat], maxclsz), sep=":"), sep=":")
    ##b <- c(rep(b[1:ncat], maxclsz), b[-(1:ncat)])
  }

  xmat <- cbind(xc, x) # note the negate sign!!!
  p <- ncol(xmat)

  offset <- model.extract(m, offset)
  if (is.null(offset)) offset <- rep(0, length(id))
  offset <- - rep(offset, rep(ncat, sum(clusz)))

  w <- model.extract(m, weights)
  if (is.null(w)) w <- rep(1, length(id))
  w <- rep(w, rep(ncat, sum(clusz)))

  CORSTRS <- c("independence", "exchangeable", "NA_ar1", "unstructured", "userdefined")
  CORSTRS.ALLOWED <- c("independence", "exchangeable", "unstructured", "userdefined")
  corstrv <- pmatch(corstr, CORSTRS.ALLOWED, -1)
  if (corstrv == -1) stop("invalid corstr.")
  corstrv <- pmatch(corstr, CORSTRS)
  corr <- list(as.integer(corstrv), maxclsz)

  if (is.null(ooffset)) ooffset <- rep(0, sum(clusz*(clusz-1)/2) * ncat^2)
  if (is.null(z)) {
    if (corstrv == 5) stop("need z matrix for userdefined corstr.") 
    else z <- genZodds(clusz, waves, corstrv, ncat)
  }

  if (length(ooffset) != sum(clusz*(clusz-1)/2) * ncat^2) stop("length(ooffset) != sum(clusz*(clusz-1)) * ncat^2 detected.")
  if (corstrv > 1 && nrow(z) != sum(clusz*(clusz-1)/2) * ncat^2) stop("nrow(z) != sum(clusz*(clusz-1)) * ncat^2 detected.")
  
  waves <- rep(waves, rep(ncat, sum(clusz)))
  if (is.null(id)) stop("ID variable not found.")

  LINKS <- c("NA_identity", "logit", "probit", "cloglog", "NA_log", "NA_inverse", "NA_fisherz", "NA_lwybc2", "NA_lwylog")
  LINKS.ALLOWED <- c("logit", "probit", "cloglog")
  mean.link.v <- pmatch(mean.link, LINKS.ALLOWED, -1)
  if (mean.link.v == -1) stop("mean.link invalid.")
  mean.link.v <- pmatch(mean.link, LINKS, -1)
  
  geestr <- list(maxwave=maxclsz,
                 mean.link=rep(mean.link.v, maxclsz),
                 variance=rep(2, maxclsz),
                 sca.link=rep(1, maxclsz),
                 cor.link=5,
                 scale.fix=as.integer(scale.fix))

  p <- ncol(xmat)
  q <- ncol(z)
  if (!is.matrix(z)) z <- as.matrix(z)

  if (is.na(b)) {
    link <- mean.link
    b <- glm.fit(xmat, Y, w, family=binomial(link))$coef
  }
  if (is.na(alpha)) alpha <- rep(0,q);
  param <- list(b, alpha, gm=rep(scale.val, 1))


  ans <- .Call("ordgee_rap", Y, xmat, offset, ooffset, w, waves, z,
               clusz, ncat, rev, geestr, corr, param, control,
               PACKAGE = "geepack")

  names(ans) <- c("beta", "alpha", "gamma", "vbeta", "valpha", "vgamma",
                  "vbeta.naiv", "valpha.naiv", "valpha.stab",
                  "vbeta.ajs", "valpha.ajs", "vgamma.ajs",
                  "vbeta.j1s", "valpha.j1s", "vgamma.j1s",
                  "vbeta.fij", "valpha.fij", "vgamma.fij",
                  "error")
  ans$xnames <- dimnames(xmat)[[2]]
  ans$zcor.names <- dimnames(z)[[2]]
  names(ans$beta) <- ans$xnames
  names(ans$alpha) <- ans$zcor.names

  ans <- c(ans, list(call=scall, clusz=clusz, control=control,
                     model=list(mean.link=mean.link,
                       variance="binomial", sca.link=NULL,
                       cor.link="log", corstr=corstr, scale.fix=scale.fix)))
  class(ans) <- "geese"
  ans
}
summary.geese <- function(object, ...) {
  mean.sum <- data.frame(estimate = object$beta,
#                         nai.se = sqrt(diag(object$vbeta.naiv)),
                         san.se = sqrt(diag(object$vbeta)),
                         ajs.se = sqrt(diag(object$vbeta.ajs)),
                         j1s.se = sqrt(diag(object$vbeta.j1s)),
                         fij.se = sqrt(diag(object$vbeta.fij)))
  mean.sum$wald <- (mean.sum$estimate / mean.sum$san.se)^2
  mean.sum$p <- 1 - pchisq(mean.sum$wald, df=1)
  rownames(mean.sum) <- object$xnames
  
  corr.sum <- data.frame(estimate = object$alpha,
#                         nai.se = sqrt(diag(object$valpha.naiv)),
                         san.se = sqrt(diag(object$valpha)),
                         ajs.se = sqrt(diag(object$valpha.ajs)),
                         j1s.se = sqrt(diag(object$valpha.j1s)),
                         fij.se = sqrt(diag(object$valpha.fij)))
  corr.sum$wald <- (corr.sum$estimate / corr.sum$san.se)^2
  corr.sum$p <- 1 - pchisq(corr.sum$wald, df=1)
  if (nrow(corr.sum) > 0) rownames(corr.sum) <- object$zcor.names

  scale.sum <- data.frame(estimate = object$gamma,
                          san.se = sqrt(diag(object$vgamma)),
                          ajs.se = sqrt(diag(object$vgamma.ajs)),
                          j1s.se = sqrt(diag(object$vgamma.j1s)),
                          fij.se = sqrt(diag(object$vgamma.fij)))
  scale.sum$wald <- (scale.sum$estimate / scale.sum$san.se)^2
  scale.sum$p <- 1 - pchisq(scale.sum$wald, df=1)
  if (!is.null(object$zsca.names)) rownames(scale.sum) <- object$zsca.names

  drop <- ifelse(c(object$control$jack, object$control$j1s, object$control$fij)== 0, TRUE, FALSE)
  if (any(drop)) {
    drop <- (3:5)[drop]
    mean.sum <- mean.sum[,-drop]
    corr.sum <- corr.sum[,-drop]
    scale.sum <- scale.sum[,-drop]
  }
  
  ans <- list(mean=mean.sum, correlation=corr.sum, scale=scale.sum,
              call=object$call, model=object$model, control=object$control,
              error=object$err, clusz=object$clusz)
  class(ans) <- "summary.geese"
  ans
}

print.geese <- function(x, digits = NULL, quote = FALSE, prefix = "", ...) {
  if(is.null(digits)) digits <- options()$digits
  else options(digits = digits)
  cat("\nCall:\n")
  dput(x$call)
  cat("\nMean Model:\n")
  cat(" Mean Link:                ", x$model$mean.link, "\n")
  cat(" Variance to Mean Relation:", x$model$variance, "\n")
  cat("\n Coefficients:\n")
  print(unclass(x$beta), digits = digits)

  if (!x$model$scale.fix) {
    cat("\nScale Model:\n")
    cat(" Scale Link:               ", x$model$sca.link, "\n")
    cat("\n Estimated Scale Parameters:\n")
    print(unclass(x$gamma), digits = digits)
  }
  else cat("\nScale is fixed.\n")
  cat("\nCorrelation Model:\n")
  cat(" Correlation Structure:    ", x$model$corstr, "\n")
  if (pmatch(x$model$corstr, "independence", 0) == 0) {
    cat(" Correlation Link:         ", x$model$cor.link, "\n")
    cat("\n Estimated Correlation Parameters:\n")
    print(unclass(x$alpha), digits = digits)
  }
  ##cat("\nNumber of observations : ", x$nobs, "\n")
  ##cat("\nMaximum cluster size   : ", x$max.id, "\n")

  cat("\nReturned Error Value:  ")
  cat(x$error, "\n")
  cat("Number of clusters:  ", length(x$clusz), "  Maximum cluster size:", max(x$clusz), "\n\n")
  invisible(x)
}

print.summary.geese <- function(x, digits = NULL,
                                quote = FALSE, prefix = "", ... ) {
  if(is.null(digits)) digits <- options()$digits
  else options(digits = digits)
  cat("\nCall:\n")
  dput(x$call)
  cat("\nMean Model:\n")
  cat(" Mean Link:                ", x$model$mean.link, "\n")
  cat(" Variance to Mean Relation:", x$model$variance, "\n")
  cat("\n Coefficients:\n")
  print(x$mean, digits = digits)

  if (x$model$scale.fix == FALSE) {
    cat("\nScale Model:\n")
    cat(" Scale Link:               ", x$model$sca.link, "\n")
    cat("\n Estimated Scale Parameters:\n")
    print(x$scale, digits = digits)
  }
  else cat("\nScale is fixed.\n")
  cat("\nCorrelation Model:\n")
  cat(" Correlation Structure:    ", x$model$corstr, "\n")
  if (pmatch(x$model$corstr, "independence", 0) == 0) {
    cat(" Correlation Link:         ", x$model$cor.link, "\n")
    cat("\n Estimated Correlation Parameters:\n")
    print(x$corr, digits = digits)
  }

  ##cat("\nNumber of observations : ", x$nobs, "\n")
  ##cat("\nMaximum cluster size   : ", x$max.id, "\n")

  cat("\nReturned Error Value:    ")
  cat(x$error, "\n")
  cat("Number of clusters:  ", length(x$clusz), "  Maximum cluster size:", max(x$clusz), "\n\n")
  invisible(x)
}
# .onLoad <- function(libname, pkgname) {
#   library.dynam("geepack", pkgname, libname)
# }

.First.lib <- function(lib, pkg) {
    library.dynam("geepack", pkg, lib)
}

# .onUnload <- function(libpath) {
#   library.dynam.unload("geepack", libpath)
# }
