glm.amcmc <- function(formula, family="binomial", weights, data = parent.frame(), amcmc,
                      starting, tuning, priors, cov.model, 
                      n.samples, sub.samples, verbose=TRUE, n.report=100, ...){
  
  ####################################################
  ##Check for unused args
  ####################################################
  formal.args <- names(formals(sys.function(sys.parent())))
  elip.args <- names(list(...))
  for(i in elip.args){
    if(! i %in% formal.args)
      warning("'",i, "' is not an argument")
  }

  ####################################################
  ##formula
  ####################################################
  if(missing(formula)){stop("error: formula must be specified")}
  
  if(class(formula) == "formula"){
    
    holder <- parseFormula(formula, data)
    Y <- holder[[1]]
    X <- as.matrix(holder[[2]])
    x.names <- holder[[3]]

  }else{
    stop("error: formula is misspecified")
  }

  p <- ncol(X)
  n <- nrow(X)
    
  ##make sure storage mode is correct
  storage.mode(Y) <- "double"
  storage.mode(X) <- "double"
  storage.mode(p) <- "integer"
  storage.mode(n) <- "integer"

  ####################################################
  ##family
  ####################################################
  if(!family %in% c("binomial","poisson"))
    stop("error: family must be binomial or poisson")

  ##default for binomial
  if(family=="binomial"){
    if(missing(weights)){weights <- rep(1, n)}
    if(length(weights) != n){stop("error: weights vector is misspecified")}
  }else{
    weights <- 0
  }
  storage.mode(weights) <- "integer"
  
  ####################################################
  ##sampling method
  ####################################################
  n.batch <- 0
  batch.length <- 0
  accept.rate <- 0
    
  if(missing(amcmc)){
    
    is.amcmc <- FALSE
    if(missing(n.samples)){stop("error: n.samples need to be specified")}
    
  }else{
    
    is.amcmc <- TRUE

    names(amcmc) <- tolower(names(amcmc))
    
    if(!"n.batch" %in% names(amcmc)){stop("error: n.batch must be specified in amcmc list")}
    n.batch <- amcmc[["n.batch"]]
    
    if(!"batch.length" %in% names(amcmc)){stop("error: batch.length must be specified in amcmc list")}
    batch.length <- amcmc[["batch.length"]]

    if(!"accept.rate" %in% names(amcmc)){
      warning("accept.rate was not specified in the amcmc list and was therefore set to the default 0.43")
      accept.rate <- 0.43
    }else{
      accept.rate <- amcmc[["accept.rate"]]
    }

    n.samples <- n.batch*batch.length
  }
  
  storage.mode(n.samples) <- "integer"
  storage.mode(n.batch) <- "integer"
  storage.mode(batch.length) <- "integer"
  storage.mode(accept.rate) <- "double"
  
  ####################################################
  ##Starting values
  ####################################################
  beta.starting <- 0

  if(missing(starting)){stop("error: starting value list for the parameters must be specified")}
  
  names(starting) <- tolower(names(starting))   

  if("beta" %in% names(starting)){
    beta.starting <- starting[["beta"]]
    if(length(beta.starting) != p){stop(paste("error: starting values for beta must be of length ",p,sep=""))}
  }else{
    if(family=="poisson"){
      beta.starting <- coefficients(glm(Y~X-1, family="poisson"))
    }else{
      beta.starting <- coefficients(glm((Y/weights)~X-1, weights=weights, family="binomial"))
    }
  }
  
  storage.mode(beta.starting) <- "double"
 
  ####################################################
  ##Priors
  ####################################################
  beta.Norm <- 0
  beta.prior <- "flat"

  if(missing(priors)) {stop("error: prior list for the parameters must be specified")}
    
  names(priors) <- tolower(names(priors))

  if("beta.normal" %in% names(priors)){
    beta.Norm <- priors[["beta.normal"]]
    if(!is.list(beta.Norm) || length(beta.Norm) != 2){stop("error: beta.Norm must be a list of length 2")}
    if(length(beta.Norm[[1]]) != p ){stop(paste("error: beta.Norm[[1]] must be a vector of length, ",p, " with elements corresponding to betas' mean",sep=""))}
    if(length(beta.Norm[[2]]) != p ){stop(paste("error: beta.Norm[[2]] must be a vector of length, ",p, " with elements corresponding to betas' sd",sep=""))}
    beta.prior <- "normal"
  }

  ####################################################
  ##Tuning values
  ####################################################
  beta.tuning <- 0

  if(!missing(tuning)){
    
    names(tuning) <- tolower(names(tuning))
    
    if(!"beta" %in% names(tuning)){stop("error: beta must be specified in tuning value list")}
    beta.tuning <- tuning[["beta"]]
    
    if(is.matrix(beta.tuning)){
      if(nrow(beta.tuning) != p || ncol(beta.tuning) != p)
        stop(paste("error: if beta tuning is a matrix, it must be of dimension ",p,sep=""))
      
      if(is.amcmc)
        beta.tuning <- diag(beta.tuning)
      
    }else if(is.vector(beta.tuning)){
      if(length(beta.tuning) != p)
        stop(paste("error: if beta tuning is a vector, it must be of length ",p,sep=""))
      
      if(!is.amcmc && p > 1)
        beta.tuning <- diag(beta.tuning)
      
    }else{
      stop("error: beta tuning is misspecified")
    }
    
  }else{##no tuning provided
    
    if(!is.amcmc){
      stop("error: tuning value list must be specified")
    }
    
    beta.tuning <- rep(0.01,p)
    
  }
  
  storage.mode(beta.tuning) <- "double"

  
  ####################################################
  ##Other stuff
  #################################################### 
  if(missing(sub.samples)){sub.samples <- c(1, n.samples, 1)}
  if(length(sub.samples) != 3 || any(sub.samples > n.samples) ){stop("error: sub.samples misspecified")}
  
  storage.mode(n.report) <- "integer"
  storage.mode(verbose) <- "integer"

  ####################################################
  ##Pack it up and off it goes
  ####################################################

  out <- .Call("glm_amcmc", Y, X, p, n, family, weights,
               beta.prior, beta.Norm, beta.starting, beta.tuning,
               n.batch, batch.length, accept.rate, verbose, n.report)

  out$weights <- weights
  out$family <- family
  out$Y <- Y
  out$X <- X
  out$n <- n
  out$p <- p
  out$verbose <- verbose
  out$sub.samples <- sub.samples
 
  ##subsample 
  out$p.samples <- mcmc(t(out$p.samples[,seq(sub.samples[1], sub.samples[2], by=as.integer(sub.samples[3]))]))
  out$n.samples <- nrow(out$p.samples)##get adjusted n.samples
  colnames(out$p.samples) <-  x.names

  ####################################################
  ##DIC
  ####################################################
  n.samples <- out$n.samples
  status <- 0
  
  d <- rep(0, n.samples)
  DIC <- matrix(0,4,1)
  
  beta <- as.matrix(out$p.samples)
  beta.mu <- colMeans(beta)
  
  for(s in 1:n.samples){
    
    if(family == "poisson"){
      d[s] <- -2*sum(-exp(X%*%beta[s,])+Y*(X%*%beta[s,]))
    }else if(family == "binomial"){
      pp <- 1/(1+exp(-X%*%beta[s,]))
      d[s] <- -2*sum(Y*log(pp)+(weights-Y)*log(1-pp))
    }else{
      stop("error: family is misspecified")
    }
  }
  
  d.bar <- mean(d)
  
  if(family == "poisson"){
    d.bar.omega <- -2*sum(-exp(X%*%beta.mu)+Y*(X%*%beta.mu))
  }else if(family == "binomial"){
    pp <- 1/(1+exp(-X%*%beta.mu))
    d.bar.omega <- -2*sum(Y*log(pp)+(weights-Y)*log(1-pp))
  }else{
    stop("error: family is misspecified")
  }
  
  pd <- d.bar - d.bar.omega
  dic <- d.bar + pd
  
  DIC <- matrix(0,4,1)
  
  rownames(DIC) <- c("bar.D", "D.bar.Omega", "pD", "DIC")
  colnames(DIC) <- c("value")
  DIC[1,1] <- d.bar
  DIC[2,1] <- d.bar.omega
  DIC[3,1] <- pd
  DIC[4,1] <- dic
  
  out$DIC <- DIC
  
  ##out it goes
  class(out) <- "glm"
  out  
}