R/CID-gibbs.R

In bdabbs13/CIDnetworks: Generative Models for Complex Networks with Conditionally Independent Dyadic Structure

#####################################################################################
                                        #
                                        # Gibbs Sampler collection for CID, given the collection of input terms.

                                        #source("COV-reference.R"); source("SBM-reference.R"); source("LSM-reference.R"); source("SR-reference.R"); library(Rcpp); library(mvtnorm); library(msm); sourceCpp ("../src/cid.cpp"); source("CID-basefunctions.R");

.onAttach <- function (...) {
    packageStartupMessage("CIDnetworks v 0.8.0")
}

                                        # All subclasses available to CID.

LSM <- function(...) LSMcid$new(...)
LVM <- function(...) LVMcid$new(...)
SBM <- function(...) SBMcid$new(...)
SR <- function(...) SRcid$new(...)
BETA <- function(...) BETAcid$new(...)


EdgeCOV <- function(..., cov.type="Edge") COVARIATEcid$new(..., cov.type=cov.type)
SenderCOV <- function(...) COVARIATEcid$new(..., cov.type="Sender")
ReceiverCOV <- function(...) COVARIATEcid$new(..., cov.type="Receiver")
SendRecCOV <- function(...) COVARIATEcid$new(..., cov.type="SendRec")
IdenticalCOV <- function(...) COVARIATEcid$new(..., cov.type="Identical")


MMSBM <- function(...) MMSBMcid$new(...)
HBM <- function(...) HBMcid$new(...)


unwrap.CID.Gibbs <- function (gibbs.out) list.output.to.matrices(gibbs.out)

                                        #{  elements <- lapply (1:length(gibbs.out[[1]]), function(el) if (grepl("out", class(gibbs.out[[1]][[el]]))) {    r1 <- lapply(1:length(gibbs.out[[1]][[el]]), function (el2) sapply(gibbs.out, function(it) it[[el]][[el2]]))    names(r1) <- names(gibbs.out[[1]][[el]])    r1  } else sapply(gibbs.out, function(it) it[[el]]))  names(elements) <- names(gibbs.out[[1]])  elements}



CIDnetwork <-
    setRefClass (
        "CIDnetwork",
        fields = list(
            n.nodes="numeric",
            edge.list="matrix",
            edge.list.rows="list",
            outcome="numeric",
            node.names="character",

            is.directed="logical",
            class.outcome="character",   #"ordinal" or "gaussian" -- if "binary", it will make it "ordinal"
            ordinal.count="numeric",
            ordinal.cutoffs="numeric",   #current draw. Every edge has the same cutoff value with respect to the outcome.

            int.outcome="numeric",  #intermediate outcome, which is Gaussian. If outcome is binary, this is what is truncated.

            intercept="numeric",    #current draw
            intercept.m="numeric",  #prior mean
            intercept.v="numeric",  #prior variance

            residual.variance="numeric",   #current draw -- always 1 if binary/ordinal
            residual.variance.ab="numeric",  #inverse gamma prior parameters -- length 2

            reciprocity.present="logical",
            reciprocal.match="integer",   #which edge has the other as its counterpart? Use for reciprocity.
            int.correlation="numeric",
            int.correlation.ab="numeric",  #symmetric beta prior.

            ## robit.augment="numeric",

            log.likelihood="numeric",   #current log likelihood of data

            components="list",   # updated with every draw -- contains COV, LSM, etc.
            comp.values="matrix" # after each draw, updates the value as expressed of that object. "value" is the component of the linear predictor for that edge ; each row corresponds to each edge; each column corresponds to each component.
        ),

        methods=list(
            initialize = function (

                                   edge.list,    #specs: this should be numbered from 1 to the max edge number.
                                        #        sociomatrix,

                                   edge.list.rows,
                                   n.nodes,      #this should come from the data.
                                   node.names=character(),

                                   intercept=0,
                                   intercept.m=0,
                                   intercept.v=1000000,

                                   residual.variance=1,
                                   residual.variance.ab=c(0.001, 0.001),

                                   outcome=numeric(0),
                                   generate=FALSE,

                                   is.directed,
                                   class.outcome="ordinal",
                                   ordinal.count=2,
                                   ordinal.cutoffs=sort(rexp(ordinal.count-2, rate=50)),

                                   int.correlation=0,          #For now, generate nothing.
                                   int.correlation.ab=c(1,1),  #Symmetric beta prior.
                                   include.reciprocity=FALSE,  #No reciprocity unless asked.

                                   reinit=TRUE,
                                   components=list(),
                                   verbose=2
                                   )
            {

                ## Initial checks. This is largely redundant right now, but that's OK.
                ## if (is.null(class.outcome)) class.outcome <- "ordinal"
                                        # browser()
                if (!(class.outcome %in% c("binary","gaussian","ordinal"))) stop (paste("class.outcome",class.outcome,"is not supported. Must be one of", paste(c("binary","gaussian","ordinal"), collapse=",")))

                if (int.correlation != 0) include.reciprocity <- TRUE

                if (missing(edge.list)) {
                    if (missing(n.nodes)){
                        stop ("CIDnetwork: Not detected: n.nodes, edge.list, or sociomatrix")
                    }else {

                        .self$n.nodes <<- n.nodes
                        if (!missing(is.directed))
                            .self$is.directed <<- is.directed
                        else {
                            .self$is.directed <<- include.reciprocity | (int.correlation != 0)
                            if(verbose > 0){
                                if (.self$is.directed)
                                    message ("CIDnetwork: Reciprocity was specified: assuming a DIRECTED network.")
                                else
                                    message ("CIDnetwork: Directedness was unspecified: assuming an undirected network.")
                            }
                        }

                        if (.self$is.directed) {
                            .self$edge.list <<- make.arc.list (n.nodes)
                            .self$outcome <<- outcome
                        } else {
                            .self$edge.list <<- make.edge.list (n.nodes)
                            .self$outcome <<- outcome
                        }
                    }
                } else {
                    if (missing(n.nodes))
                        .self$n.nodes <<- max(c(edge.list))
                    else
                        .self$n.nodes <<- n.nodes

                    if (missing(is.directed)) {

                        e1 <- unique(edge.list)

                        e1[e1[,1]>e1[,2],] <- e1[e1[,1]>e1[,2], 2:1]
                        e2 <- unique(e1)
                        if (nrow(e1) != nrow(e2)) {
                            if(verbose > 0)
                                message ("CIDnetwork: Reciprocal potential edges detected; assuming this is a directed network.")
                            .self$is.directed <<- TRUE
                        } else {
                            if(verbose > 0)
                                message ("CIDnetwork: No reciprocal potential edges detected; assuming this is an undirected network.")
                            .self$is.directed <<- FALSE
                        }
                    }else{
                        .self$is.directed <<- is.directed
                    }

                    .self$edge.list <<- edge.list
                    .self$outcome <<- outcome
                }


                if (!missing(edge.list.rows))
                    .self$edge.list.rows <<- edge.list.rows
                else
                    .self$edge.list.rows <<- row.list.maker(.self$edge.list)

                if (length(node.names) != .self$n.nodes){
                    if(verbose > 0)
                        message("Warning:  Length of node.names differs from n.nodes.  Using default node names")
                    .self$node.names <<- as.character(1:.self$n.nodes)
                }else{
                    .self$node.names <<- node.names
                }
                if (class.outcome == "binary") {
                    .self$class.outcome <<- "ordinal"
                    .self$ordinal.count <<- 2
                    .self$ordinal.cutoffs <<- numeric()
                } else {
                    .self$class.outcome <<- class.outcome
                                        #          if(class.outcome != "gaussian"){
                    .self$ordinal.count <<- ordinal.count
                    .self$ordinal.cutoffs <<- ordinal.cutoffs
                                        #          }
                }

                .self$intercept <<- intercept
                .self$intercept.m <<- intercept.m
                .self$intercept.v <<- intercept.v

                if (class.outcome == "gaussian") {
                    .self$residual.variance <<- residual.variance
                } else {
                    .self$residual.variance <<- 1
                }
                .self$residual.variance.ab <<- residual.variance.ab


                ##reinitialize components, just in case.
                if (length(components)>0) {
                    if (class(components) != "list") components.t <- list(components) else components.t <- components
                    if(reinit){
                        for (kk in 1:length(components.t)) {
                            components.t[[kk]]$reinitialize(n.nodes=.self$n.nodes,
                                                            edge.list=.self$edge.list,
                                                            node.names=.self$node.names)

                            if (class(components.t[[kk]]) %in% c("SBMcid", "MMSBMcid", "HBMcid")) {
                                .self$intercept.m <<- 0
                                .self$intercept.v <<- 0.000000000001
                            }
                        }
                    }
                    .self$components <<- components.t
                                        #          if(length(which(sapply(components ,class) == "COVARIATEcid")) >1){
                                        #            .self$components <- combine.covariates(components)
                                        #          }
                } else .self$components <- list()

                ##message("Component initialization complete.")



                ##Let's do reciprocity now! ACT, 5-21-14
                .self$int.correlation <<- int.correlation
                .self$int.correlation.ab <<- int.correlation.ab
                if(verbose > 1) message("Reciprocity: ",include.reciprocity)

                reciprocity.present <<- include.reciprocity
                if (reciprocity.present) {
                    ## find matches.
                    reciprocal.match <<- apply(.self$edge.list, 1, function(ee) min(which (.self$edge.list[,1] == ee[2] & .self$edge.list[,2] == ee[1])))
                } else {
                    reciprocal.match <<- rep(as.integer(NA), length(outcome))
                }
                ##int.correlation=0,          #For now, generate nothing.
                ##int.correlation.ab=c(1,1),  #Symmetric beta prior.
                ##include.reciprocity=FALSE,  #No reciprocity unless asked.



                if (generate) .self$generate()

                ## Control statements currently force .self$outcome to have already
                ## been set anyways.
                ##        else if (length(outcome) == nrow(.self$edge.list)) {
                ##          .self$outcome <<- outcome
                ##        }

                if (length(.self$outcome) != nrow(.self$edge.list))
                    stop (paste0("In CIDnetwork$initialize: Outcome variable length: ",length(.self$outcome),"; Edges: ", nrow(.self$edge.list),". Did you forget to add a required variable?"))

                ## Check for empty categories, currently disallowed.
                if (class.outcome == "ordinal") {
                    counts <- tabulate (.self$outcome + 1, .self$ordinal.count)
                    if (any(counts == 0))
                        stop ("The following ordinal categories have no outcomes: ", paste(which(counts == 0)-1, collapse=" "), "; required non-zero counts for ", paste (0:(.self$ordinal.count-1), collapse=" "))
                }

                ##message("Generation complete.")
                update.intermediate.outcome()
                ##message("CID Initialization complete.")

            },

            reinitialize = function (n.nodes=NULL, edge.list=NULL, node.names=NULL) {
                if (!is.null(n.nodes)) n.nodes <<- n.nodes
                if (!is.null(edge.list)) {
                    edge.list <<- edge.list
                    edge.list.rows <<- row.list.maker(edge.list)
                }

                if (!is.null(node.names)) {
                    if (length(node.names) == .self$n.nodes) node.names <<- node.names
                } else node.names <<- as.character(1:.self$n.nodes)


                if (length(components)>0) for (kk in 1:length(components)) {
                                              components[[kk]]$reinitialize (.self$n.nodes, .self$edge.list, .self$node.names)
                                              if (class(components[[kk]]) %in% c("SBMcid", "MMSBMcid", "HBMcid")) {
                                                  intercept.m <<- 0
                                                  intercept.v <<- 0.000000000001
                                              }
                                          }


                if (reciprocity.present) {
                    ## find matches.
                    reciprocal.match <<- apply(edge.list, 1, function(ee) min(which (edge.list[,1] == ee[2] & edge.list[,2] == ee[1])))
                } else {
                    reciprocal.match <<- rep(as.integer(NA), length(outcome))
                }

            },

            pieces = function (include.name=FALSE) {
                if (length(components)>0) {
                    c(list(intercept=intercept,
                           residual.variance=residual.variance,
                           log.likelihood=log.likelihood,
                           ordinal.cutoffs=ordinal.cutoffs,
                           int.correlation=int.correlation),
                      lapply(components, function(cc) cc$pieces(include.name)))
                } else {
                    list(intercept=intercept,
                         residual.variance=residual.variance,
                         log.likelihood=log.likelihood,
                         ordinal.cutoffs=ordinal.cutoffs,
                         int.correlation=int.correlation)
                }
            },

            show = function () {
                message("CIDnetwork object properties:")

                message(paste("class.outcome:", class.outcome))
                if (class.outcome == "ordinal") message(paste("Ordinal groups:", ordinal.count))
                message(paste("Nodes:", n.nodes))
                if (is.directed) message(paste("Potential Directed Arcs:", nrow(edge.list))) else message(paste("Potential Undirected Edges:", nrow(edge.list)))


                message("Intercept: ", intercept)
                message("Variance: ", residual.variance)
                if (length(components)>0) for (kk in 1:length(components)) {
                                              message(class(components[[kk]])); components[[kk]]$show()
                                          }

            },
                                        #      combine.covariates <- function(components.arg){
                                        #        cov.ind <- which(sapply(components,class) == "COVARIATEcid")
                                        #
                                        #    },
            plot = function (coefs=coef.cov, names=1:length(coefs), sd=NULL, interval=NULL, ...) {
                if (length(components) > 0) for (cc in 1:length(components)) components[[cc]]$plot()
            },
            plot.network = function (color=outcome, ...) {
                image.netplot(edge.list, color, node.labels=node.names, ...)
            },

            value = function (redo=FALSE) {
                if (redo) update.comp.values()
                rowSums(comp.values) + intercept
            },

            value.ext = function(parameters=pieces(), edges=1:nrow(edge.list)) {
                if (length(components) > 0){
                    ncc <- non.comp.count()
                    # comp.seq <- seq(non.comp.count()+1,non.comp.count() + length(components))
                    # par_components <- parameters[comp.seq]
                    comp.values.here <- parameters$intercept
                    for(cc in 1:length(components)){
                        comp.values.here <- cbind(comp.values.here,
                                                  components[[cc]]$value.ext(parameters[[ncc+cc]]))
                    }
                                        #comp.values.here <- sapply (par_components, function(cc) cc$value())
                }else{
                    comp.values.here <- matrix(parameters$intercept,nrow=nrow(edge.list))
                }
                # (rowSums(comp.values.here) + parameters$intercept)[edges]
                rowSums(comp.values.here)[edges]
            },

            summary = function () show(),


            generate = function () {

                                        # Built 5-21-14. Tested: ----
                if (reciprocity.present) {
                    int.deviation <- matrix(rnorm(2*nrow(edge.list), 0, sqrt(residual.variance)), ncol=2)
                    int.deviation[,2] <- int.correlation*int.deviation[,1] + sqrt(1-int.correlation^2)*int.deviation[,2]

                    swaps <- which(edge.list[,1]>edge.list[,2] & !is.na(reciprocal.match))

                    int.deviation[swaps,1] <- int.deviation[reciprocal.match[swaps],2]
                    int.outcome <<- value(redo=TRUE) + int.deviation[,1]
                } else {
                    int.outcome <<- rnorm(nrow(edge.list),
                                          value(redo=TRUE),
                                          sqrt(residual.variance))
                }

                ## Generate user-observed outcome from intermediate outcome.
                if (class.outcome == "ordinal") {
                    temp.out <- 0*int.outcome
                    ordinal.steps <- c(0, ordinal.cutoffs)
                    for (ii in 1:length(ordinal.steps)) temp.out <- temp.out + 1*(int.outcome > ordinal.steps[ii])
                    outcome <<- temp.out
                }

                if (class.outcome == "gaussian") outcome <<- int.outcome

            },


            rem.values = function(kk) {if (kk>0) value() - comp.values[,kk] else value() - intercept},

            update.comp.values = function () {
                if (length(components)>0){
                    comp.values <<- sapply (components, function(cc) cc$value())
                }else{
                    comp.values <<- matrix(0, nrow=nrow(edge.list))
                }
            },



            ## Update Z_ij.
            update.intermediate.outcome = function () {
                value.hold <- value(redo=TRUE)
                if(class.outcome == "gaussian"){
                    int.outcome <<- outcome
                }
                if(any(is.na(outcome))){
                    io.temp <- int.outcome
                    io.temp[is.na(outcome)] <- rnorm(sum(is.na(outcome)),
                                                     value.hold[is.na(outcome)],1)
                    int.outcome <<- io.temp

                }
                if (class.outcome == "ordinal") {
                    ##first, update the values themselves.
                    io.temp <- int.outcome

                    ##matched.value <- value.hold[reciprocal.match]
                    lohi <- edge.list[,1] < edge.list[,2]

                    breaker.lower <- c(-Inf, 0, ordinal.cutoffs)
                    breaker.upper <- c(0, ordinal.cutoffs, Inf)

                    breakers.lower <- breaker.lower[outcome + 1]
                    breakers.upper <- breaker.upper[outcome + 1]


                    ##are any going to be trouble?
                    p.gen <- rep(0.5, length(outcome))
                    for (ii in 1:ordinal.count) {
                        p.gen[!is.na(outcome) & outcome == ii-1] <-
                            pnorm(breaker.upper[ii], value.hold[!is.na(outcome) & outcome == ii-1], 1) - pnorm(breaker.lower[ii], value.hold[!is.na(outcome) & outcome == ii-1], 1)
                    }
                    p.gen[is.na(outcome)] <- 0

                    ##for (ii in 1:ordinal.count) {
                    ## unmatched edges.
                    condition <- p.gen > 1e-10 & is.na(reciprocal.match)
                    io.temp[condition] <-
                        rtnorm(sum(condition),
                               value.hold[condition], 1,
                               lower=breakers.lower[condition],
                               upper=breakers.upper[condition])

                    ##matched edges, (lower, higher).
                    condition <- p.gen > 1e-10 & !is.na(reciprocal.match) & lohi
                    io.temp[condition] <-
                        rtnorm(sum(condition),
                               value.hold[condition] + int.correlation*(io.temp[reciprocal.match[condition]] - value.hold[reciprocal.match[condition]]),
                               1 - int.correlation^2,
                               lower=breakers.lower[condition],
                               upper=breakers.upper[condition])

                                        #matched edges, (higher, lower).
                    condition <- p.gen > 1e-10 & !is.na(reciprocal.match) & !lohi
                    io.temp[condition] <-
                        rtnorm(sum(condition),
                               value.hold[condition] + int.correlation*(io.temp[reciprocal.match[condition]] - value.hold[reciprocal.match[condition]]),
                               1 - int.correlation^2,
                               lower=breakers.lower[condition],
                               upper=breakers.upper[condition])

                    p.gen[is.na(outcome)] <- 0.5
                    io.temp[p.gen <= 1e-10] <- (breakers.lower[p.gen <= 1e-10]+breakers.upper[p.gen <= 1e-10])/2
                    io.temp[p.gen <= 1e-10 & outcome == 0 ] <- (breakers.upper[outcome == 0 & p.gen <= 1e-10])
                    io.temp[(p.gen <= 1e-10) & (outcome == ordinal.count - 1)] <- (breakers.lower[outcome == ordinal.count - 1 & p.gen <= 1e-10])

                    int.outcome <<- io.temp

                    ##now, change the cutoff values, which lie between the Z values for each one. Assume a flat prior for now.
                    if (length(ordinal.cutoffs)>0) for (kk in 1:length(ordinal.cutoffs)) {
                                                       effective.range <- c(max(int.outcome[outcome <= kk]), min(int.outcome[outcome >= kk+1]))
                                                       if (is.na(effective.range[1])) effective.range[1] <- 0
                                                       if (is.na(effective.range[2])) effective.range[1] <- 10000
                                                       ordinal.cutoffs[kk] <<- runif(1, effective.range[1], effective.range[2])
                                                   }

                }
            },



            ## Simple slice sampler for autocorrelation term \rho. Pretty optimized! Act, 6-3-14
            int.correlation.prior = function (pp = int.correlation) {
                dbeta((1+pp)/2,  int.correlation.ab[1], int.correlation.ab[2], log=TRUE) - log(2)
            },

                                        #Note: this is being drawn from the main likelihood, not the intermediate likelihood, because it's more efficient for the whole chain.
            draw.int.correlation = function () {

                these.pieces <- pieces()
                current.value <- int.correlation

                                        #First: draw uniform at current value and establish the vertical cut.
                cut.1 <- log(runif(1)) +
                    log.likelihood.by.value(value.ext(these.pieces), these.pieces) + #, use.intermediate=TRUE) +  BD 12/31
                    int.correlation.prior(current.value)
                limits <- c(-0.9999,0.9999)

                emergency.count <- 0
                repeat {
                    prop.value <- runif(1, limits[1], limits[2])   #Draw proposal value between the bounds.
                    these.pieces$int.correlation <- prop.value
                    if (log.likelihood.by.value(value.ext(these.pieces), these.pieces) +  #, use.intermediate=TRUE) + BD
                        int.correlation.prior(prop.value) > cut.1) {   #is the current density above the threshold?
                        int.correlation <<- prop.value; break    #keep it and save.
                    } else limits[1 + 1*(prop.value > current.value)] <- prop.value    #Trim down.

                    emergency.count <- emergency.count+1; if (emergency.count > 1000) stop ("In draw.int.correlation, way too many narrows-down were made. This is either at a peak or it's broken.")
                }

            },

                                        #What does the grid say?
            test.int.correlation = function (rho=seq(-19,19)/20) {
                these.pieces <- pieces()
                sapply (rho, function (rr) {
                    these.pieces$int.correlation <- rr
                    log.likelihood.by.value(value.ext(these.pieces), these.pieces) + #, use.intermediate=TRUE) + BD 12/31
                        int.correlation.prior(rr)
                })
            },




            draw.intercept = function (verbose=2) {

                outcomeresid <- int.outcome - rem.values(0);
                varpiece <- solve(nrow(edge.list)/residual.variance + 1/intercept.v)
                meanpiece <- varpiece*(sum(outcomeresid)/residual.variance + intercept.m/intercept.v)
                if (verbose > 2) message ("Intercept ",meanpiece," ",sqrt(varpiece))
                intercept <<- rnorm(1, meanpiece, sqrt(varpiece))

            },

            ## Note: getting this to work for 2D pmvnorm. ACT, 5-23-14
            log.likelihood.by.value = function (value.this=value(),
                                                pieces.this=pieces(),
                                                sumup=TRUE,
                                                use.intermediate=FALSE)
            {
                cm <- function(pp) matrix(c(1,pp,pp,1), nrow=2)
                output <- 0*int.outcome
                if (class.outcome == "gaussian" | use.intermediate) {
                    outcomeresid <- int.outcome - value.this
                    if (reciprocity.present) {
                        picks <- which(is.na(reciprocal.match))
                        output[picks] <- dnorm(outcomeresid[picks], 0, sqrt(pieces.this$residual.variance), log=TRUE)

                        picks2 <- which (!is.na(reciprocal.match) & edge.list[,1] < edge.list[,2])
                        output[picks2] <- dmvnorm(cbind(outcomeresid[picks2], outcomeresid[reciprocal.match[picks2]]),
                                                  rep(0, 2),
                                                  pieces.this$residual.variance*cm(pieces.this$int.correlation),
                                                  log=TRUE)
                    } else output <- dnorm(outcomeresid, 0, sqrt(pieces.this$residual.variance), log=TRUE)

                    if (sumup) output <- sum(output[!is.na(outcome)])
                } else if (class.outcome == "ordinal") {

                    breaker.lower <- c(-Inf, 0, pieces.this$ordinal.cutoffs)
                    breaker.upper <- c(0, pieces.this$ordinal.cutoffs, Inf)

                    ##  Getting ranges for non-NA outcomes
                    breakers.lower <- breaker.lower[outcome[!is.na(outcome)] + 1]
                    breakers.upper <- breaker.upper[outcome[!is.na(outcome)] + 1]
                    value.this <- value.this[!is.na(outcome)]

                    ##  BD Added 12/31
                    reciprocal.match.na <- reciprocal.match[!is.na(outcome)]

                                        # first: unmatched edges.
                    if (reciprocity.present) {
                        picks <- which(is.na(reciprocal.match.na))
                        output[picks] <- log(
                            pnorm(breakers.upper[picks],
                                  value.this[picks],
                                  sqrt(pieces.this$residual.variance)) -
                            pnorm(breakers.lower[picks],
                                  value.this[picks],
                                  sqrt(pieces.this$residual.variance)))

                        ## now, matched edges.  BD added !is.na(outcome)
                        picks2 <- which (!is.na(reciprocal.match.na) & edge.list[!is.na(outcome),1] < edge.list[!is.na(outcome),2])

                        ## ACT was working on this. See basefunctions line 51
                        output[picks2] <- my.pmvnorm(lower=matrix(breakers.lower[c(picks2, reciprocal.match.na[picks2])], ncol=2),
                                                     upper=matrix(breakers.upper[c(picks2, reciprocal.match.na[picks2])], ncol=2),
                                                     meanval=matrix(value.this[c(picks2, reciprocal.match.na[picks2])], ncol=2),
                                                     sigma=pieces.this$residual.variance,
                                                     rho=pieces.this$int.correlation)

                                        #this.output <- sapply(picks2, function (pp) {
                                        #  pmvnorm (breakers.lower[c(pp, reciprocal.match.na[pp])],
                                        #          breakers.upper[c(pp, reciprocal.match.na[pp])],
                                        #           mean=value.this[c(pp, reciprocal.match.na[pp])],
                                        #           sigma=pieces.this$residual.variance*matrix(c(1, pieces.this$int.correlation, pieces.this$int.correlation, 1), nrow=2))})
                                        #output[picks2] <- log(this.output)

                    } else {

                        output <- log(
                            pnorm(breakers.upper,
                                  value.this,
                                  sqrt(pieces.this$residual.variance)) -
                            pnorm(breakers.lower,
                                  value.this,
                                  sqrt(pieces.this$residual.variance)))
                    }
                    if(sumup) output <- sum(output)
                }

                return(output)
            },

            update.log.likelihood = function () {
                log.likelihood <<- log.likelihood.by.value ()
            },

            draw.variance = function (verbose=2) {
                outcomeresid <- int.outcome - value();

                if (verbose > 2) message ("Variance: ",nrow(edge.list), " ", sum(outcomeresid^2))
                residual.variance <<-
                    1/rgamma(1,
                             residual.variance.ab[1] + nrow(edge.list)/2,
                             residual.variance.ab[2] + sum(outcomeresid^2)/2)

            },

            draw = function (verbose=FALSE) {

                if (reciprocity.present) draw.int.correlation()
                if (class.outcome != "gaussian") update.intermediate.outcome()

                if (length(components)>0) for (kk in 1:length(components)) {
                                              update.comp.values()
                                              components[[kk]]$outcome <<- .self$int.outcome - rem.values(kk)
                                              components[[kk]]$residual.variance <<- residual.variance
                                              components[[kk]]$draw()

                                              if (exists("shift", components[[kk]])) {   ### - where does this currently arise?
                                                  intercept <<- intercept + components[[kk]]$shift
                                                  components[[kk]]$shift <<- 0
                                              }
                                          }

                                        #variance and intercept.
                update.comp.values()
                draw.intercept(verbose)

                if (class.outcome == "gaussian") {
                    update.comp.values()
                    draw.variance(verbose)
                    update.comp.values()
                }

                ## update correlation!!! ACT, 5-23-14

                update.log.likelihood()

            },

            random.start = function () {
                intercept <<- rnorm (1, 0, 1)
                                        #        draw.intercept()
                if (length(components)>0) for (kk in 1:length(components)) components[[kk]]$random.start()
                if (class.outcome == "gaussian") draw.variance()
                if (reciprocity.present) draw.int.correlation()

                update.log.likelihood()

            },

            gibbs.full = function (
                                   report.interval=100, draws=100,
                                   burnin=-1, thin=10,   #auto-cut.
                                   auto.burn.cut=TRUE,
                                   auto.burn.count=200,
                                   make.random.start=TRUE,
                                   auto.converge=FALSE,
                                   extend.max = 10,
                                   extend.count = 100,
                                   verbose = 2) {

                                        #if (auto.burn.cut) {thin <- 10; burnin <- 0}
                if (thin < 1) stop ("thin must be an integer greater than zero.")
                if (burnin >= 0) {
                    if(verbose > 1)
                        message ("Overriding auto burn-in and cut. Initially discarding ",burnin," iterations.")
                    auto.burn.cut <- FALSE
                } else if (burnin < 0 & !auto.burn.cut) {
                    if(verbose > 0)
                        message ("Negative burn-in inputted -- defaulting to auto-burn-in.")
                    auto.burn.cut <- TRUE
                }


                out <- list()
                if (make.random.start) random.start()

                if (auto.burn.cut) {
                    if(verbose > 1)
                        message ("CID Auto-Burning In")
                    repeat {
                        loglikes <- rep(NA, auto.burn.count)
                        theseruns <- 1:auto.burn.count
                        time.one <- proc.time()[3]
                        for (kk in theseruns) {
                            draw()
                            loglikes[kk] <- log.likelihood
                        }
                        time.two <- proc.time()[3] - time.one

                                        #  obj1 <- lm(loglikes ~ theseruns)
                                        #  slopes <- summary(obj1)$coef[2, c(1,4)]
                                        #  if (slopes[1]<0 | slopes[2] > 0.05) break else message ("CID Still Auto-Burning In")
                                        #  if (any(is.na(loglikes))) {print (loglikes)} else {
                        slopes <- sum(loglikes*(theseruns - mean(theseruns)))
                        if (is.na(slopes) | slopes >= 0){
                            if(verbose > 1) message ("CID Still Auto-Burning In")
                        }else ##if (slopes<0)
                            break
                        ##else message ("CID Still Auto-Burning In")
                        ##}
                    }

                    if(verbose > 0)
                        message ("CID Auto-Burned In. Estimated Seconds Remaining: ",
                                 round(time.two/auto.burn.count*thin*draws))
                    burnin <- 0
                }


                if(burnin > 0) time.one <- proc.time()[3]
                for (kk in 1:(draws*thin+burnin)) {
                    draw();
                    index <- (kk-burnin)/thin
                    if (kk > burnin & round(index)==index) {

                        out[[index]] <- pieces(include.name=TRUE)
                        if (report.interval > 0 & index %% report.interval == 0) {
                            if(verbose > 1) message("CID ",index)
                        }
                    } else if (round(index)==index) {
                        if (report.interval > 0) if (index %% report.interval == 0){
                                                     if(index != 0){
                                                         if(verbose > 1) message("CID burnin ",index)
                                                     }else{
                                                         time.two <- proc.time()[3] - time.one
                                                         if(verbose > 1) message ("CID Burned In. Estimated Seconds Remaining: ",
                                                                                  round((time.two/burnin)*thin*draws))
                                                     }
                                                 }
                    }
                }

                if(auto.converge){
                    extend.iter <- 0
                    extend.count <- min(extend.count,draws)

                    while(!convergence.test(out) & extend.iter < extend.max) {
                        if(verbose > 0)
                            message("CID Convergence not detected.  Extending chain")
                        extend.iter <- extend.iter + 1

                        ##  Shifting Chain
                        if(draws > extend.count){
                            for(ii in 1:(draws - extend.count)){
                                out[[ii]] <- out[[ii + extend.count]]
                            }
                        }

                        ##  Getting new draws
                        for(kk in 1:(extend.count*thin)){
                            draw();
                            index <- draws - extend.count + kk/thin
                            if(round(index) == index){
                                out[[index]] <- pieces(include.name=TRUE)
                            }
                        }
                    }

                    if(!convergence.test(out)){
                        if(verbose > 0)
                            message("CID Convergence not detected after maximum extensions")
                    }else{
                        if(verbose > 0)
                            message("CID Converged")
                    }
                }

                return(out)
            },

            convergence.test = function(gibbs.out){
                log.lik <- unlist(gibbs.switcheroo(gibbs.out)$log.lik)
                draws <- length(log.lik)
                chain.1 <- log.lik[1:(draws/10)]
                chain.2 <- log.lik[(draws/2 + 1):draws]
                sd.est <- sqrt(var(chain.1)/length(chain.1) +
                               var(chain.2)/length(chain.2))

                return((mean(chain.2) - mean(chain.1)) <= qnorm(0.999) * sd.est)
            },

                                        #Removes magic number later. This should be the length of (intercept, residual, loglik, cutoffs, int.correlation).
            non.comp.count = function () 5,


            gibbs.value = function (gibbs.out){
                sapply(gibbs.out, function(gg) return(value.ext (gg)))
            },


            gibbs.mean = function(gibbs.out){

                switched <- gibbs.switcheroo(gibbs.out)
                intercept.mean <- mean(unlist(switched$intercept))
                components.mean <- NULL

                if(length(components) > 0){
                    for(cc in 1:length(components)){
                        components.mean <- c(components.mean,components[[cc]]$gibbs.mean(switched[[cc+non.comp.count()]]))
                    }
                }
                if(ordinal.count > 2){
                    ordinal.cuts.mean <- sort(sapply(switched$ordinal.cuts,mean))
                }else{
                    ordinal.cuts.mean <- ordinal.cutoffs
                }

                mean.object <- CID(edge.list,outcome,
                                   intercept=intercept.mean,
                                   components=components.mean,
                                   intercept.m=intercept.m,intercept.v=intercept.v,
                                   residual.variance=residual.variance,
                                   residual.variance.ab=residual.variance.ab,
                                        #ordinal.count=ordinal.count,
                                   class.outcome=class.outcome,
                                   ordinal.cutoffs=ordinal.cuts.mean,
                                   int.correlation=int.correlation,
                                   int.correlation.ab=int.correlation.ab,
                                   include.reciprocity=reciprocity.present,
                                   reinit=FALSE,
                                   verbose=0)
                mean.object$update.log.likelihood()
                mean.object$update.comp.values()
                return(mean.object)


            },


            gibbs.switcheroo = function (gibbs.out) {   #returns the chain for each component/parameter.
                out <- lapply(1:length(gibbs.out[[1]]), function(el)
                    lapply(1:length(gibbs.out), function(el2) gibbs.out[[el2]][[el]]))
                out.names <- c("intercept","res.variance","log.lik",
                               "ordinal.cuts","reciprocity")
                s.total <- length(out)
                if(s.total > non.comp.count()){
                    sapply(out,FUN=function(x)return(class(x[[1]])))
                    comp.classes <- sapply(out,FUN=function(x) return(class(x[[1]])))
                    out.names <- c(out.names,
                                   comp.classes[(non.comp.count()+1):s.total])
                }
                names(out) <- out.names
                return(out)
            },




                                        #      gibbs.summary = function (gibbs.out) {
                                        #        switched <- gibbs.switcheroo (gibbs.out)
                                        #        s.sum <- function (int1) c(min=min(int1), max=max(int1), mean=mean(int1), sd=sd(int1), quantile(int1, c(0.025, 0.975)))
                                        #        out <- list()

                                        #        out$intercept <- s.sum(unlist(switched[[1]]))
                                        #        message ("Intercept:"); print(out$intercept)

                                        #        out$residual.variance <- s.sum(unlist(switched[[2]]))
                                        #        if (out$residual.variance[4] != 0) {  #SD
                                        #          message ("Residual variance:")
                                        #          print (out$residual.variance)
                                        #        }

                                        #        out$log.likelihood <- s.sum(unlist(switched[[3]]))
                                        #        message ("Log likelihood:"); print(out$log.likelihood)

                                        #        if (class.outcome == "ordinal") if (ordinal.count > 2) {
                                        #          o1 <- matrix(unlist(switched[[4]]), nrow=ordinal.count-2)
                                        #          out$ordinal.cutoffs <- t(apply(o1, 1, s.sum))
                                        #          rownames(out$ordinal.cutoffs) <- paste("Cutoff", 1:(ordinal.count-2), 2:(ordinal.count-1), sep="-")
                                        #          message ("Ordinal cutoffs:")
                                        #          print (out$ordinal.cutoffs)
                                        #        }

                                        #        if (length(components) > 0) for (cc in 1:length(components)) {
                                        #          message ("Component ", class(components[[cc]]),":")
                                        #          out[[cc+non.comp.count()]] <- components[[cc]]$print.gibbs.summary(switched[[cc+non.comp.count()]])
                                        #          names(out)[cc+non.comp.count()] <- class(components[[cc]])
                                        #        }

                                        #        out <- structure(out,class="summary.CID.Gibbs")
                                        #        return(invisible(out))
                                        #      },

            gibbs.summary = function (gibbs.out) {
                switched <- gibbs.switcheroo (gibbs.out$results)
                s.sum <- function (int1) {
                    c(min=round(min(int1),3), max=round(max(int1),3),
                      estimated.mean= round(mean(int1),3),
                      estimated.sd=round(sd(int1),3),
                      round(quantile(int1, c(0.025, 0.975)),3))
                }
                out <- list()

                out$intercept <- s.sum(unlist(switched[[1]]))

                out$residual.variance <- s.sum(unlist(switched[[2]]))

                out$log.likelihood <- s.sum(unlist(switched[[3]]))

                if (class.outcome == "ordinal") if (ordinal.count > 2) {
                                                    o1 <- matrix(unlist(switched[[4]]), nrow=ordinal.count-2)
                                                    out$ordinal.cutoffs <- t(apply(o1, 1, s.sum))
                                                    rownames(out$ordinal.cutoffs) <- paste("Cutoff", 1:(ordinal.count-2), 2:(ordinal.count-1), sep="-")
                                                }else{
                                                    out$ordinal.cutoffs <- NULL
                                                }

                if (length(components) > 0) for (cc in 1:length(components)) {
                                                ncc <- non.comp.count()
                                                out[[cc+ncc]] <- components[[cc]]$gibbs.summary(switched[[cc+ncc]])
                                                names(out)[cc+non.comp.count()] <- class(components[[cc]])
                                            }

                out$CID.object <- gibbs.out$CID.object
                out <- structure(out,class="summary.CID.Gibbs")
                return(out)
            },

                                        #      print.gibbs.summary = function (gibbs.out) {
                                        #        switched <- gibbs.switcheroo (gibbs.out)
                                        #        s.sum <- function (int1) {c(min=round(min(int1),3), max=round(max(int1),3), estimated.mean = round(mean(int1)), estimated.sd=round(sd(int1),3), round(quantile(int1, c(0.025, 0.975)),3))}
                                        #        out <- list()
                                        #
                                        #        out$intercept <- s.sum(unlist(switched[[1]]))
                                        #        message ("Intercept:"); print(out$intercept)
                                        #
                                        #        out$residual.variance <- s.sum(unlist(switched[[2]]))
                                        #        if (out$residual.variance[4] != 0) {  #SD
                                        #          message ("Residual variance:")
                                        #          print (out$residual.variance)
                                        #        }
                                        #
                                        #        out$log.likelihood <- s.sum(unlist(switched[[3]]))
                                        #        message ("Log likelihood:"); print(out$log.likelihood)
                                        #
                                        #        if (class.outcome == "ordinal") if (ordinal.count > 2) {
                                        #          o1 <- matrix(unlist(switched[[4]]), nrow=ordinal.count-2)
                                        #          out$ordinal.cutoffs <- t(apply(o1, 1, s.sum))
                                        #          rownames(out$ordinal.cutoffs) <- paste("Cutoff", 1:(ordinal.count-2#), 2:(ordinal.count-1), sep="-")
                                        #          message ("Ordinal cutoffs:")
                                        #          print (out$ordinal.cutoffs)
                                        #        }
                                        #
                                        #        if (length(components) > 0) for (cc in 1:length(components)) {
                                        #          message ("Component ", class(components[[cc]]),":")
                                        #          out[[cc+non.comp.count()]] <- components[[cc]]$print.gibbs.summary(#switched[[cc+non.comp.count()]])
                                        #          names(out)[cc+non.comp.count()] <- class(components[[cc]])
                                        #        }
                                        #
                                        #        out <- structure(out,class="summary.CID.Gibbs")
                                        #        return(invisible(out))
                                        #      },


            print.gibbs.summary = function(gibbs.sum){

                message ("Intercept:"); print(gibbs.sum$intercept[-(1:2)])
                if (gibbs.sum$residual.variance[4] != 0) {  #SD
                    message ("Residual variance:")
                    print (gibbs.sum$residual.variance)
                }

                message ("Log likelihood:"); print(gibbs.sum$log.likelihood[-(1:2)])

                if(!is.null(gibbs.sum$residual.variance[-(1:2)])){
                    if(gibbs.sum$residual.variance["estimated.sd"] != 0){
                        message ("Residual variance:")
                        print (gibbs.sum$residual.variance)
                    }
                }

                if (class.outcome == "ordinal") if (ordinal.count > 2) {
                                                    message ("Ordinal cutoffs:")
                                                    print (gibbs.sum$ordinal.cutoffs)
                                                }

                if(!is.null(gibbs.sum$ordinal.cutoffs)){
                    message ("Ordinal cutoffs:")
                    print (gibbs.sum$ordinal.cutoffs)
                }

                if (length(components) > 0){
                    cov.count <- 0
                    for (cc in 1:length(components)) {
                        ##		message ("Component ", class(components[[cc]]),":")
                        ix <- which(names(gibbs.sum) == class(components[[cc]]))
                        if(length(ix) > 1){
                            cov.count <- cov.count + 1
                            ix <- ix[cov.count]
                        }
                        for(ii in 1:length(ix)){
                            components[[cc]]$print.gibbs.summary(gibbs.sum[[ix]])
                        }

                                        #		print(round(gibbs.sum$class(components[[cc]])),3)
                    }
                }
                return()#invisible(gibbs.sum))
            },

            gibbs.plot = function (gibbs.out, DIC=NULL, which.plots=1:(length(components)+5), auto.layout=TRUE) {

                switched <- gibbs.switcheroo (gibbs.out)

                if (auto.layout) {
                    if (intercept.v < 0.0001) which.plots <- which.plots[which.plots != 1]
                    if (class.outcome != "gaussian") which.plots <- which.plots[which.plots != 2]
                    if (ordinal.count == 2) which.plots <- which.plots[which.plots != 4]
                    if (!reciprocity.present) which.plots <- which.plots[which.plots != 5]
                    cols <- ceiling(length(which.plots)/2)
                    par(mfrow=c(2, cols))
                }

                ## 1: Grand Intercept
                if (1 %in% which.plots & intercept.v >= 0.0001) plot.default (unlist(switched[[1]]),
                                                                              main="Grand Intercept",
                                                                              xlab="Iteration",
                                                                              ylab="Intercept",
                                                                              type="l")

                ## 2: Residual variance. Skipped if not gaussian.
                if (2 %in% which.plots & class.outcome=="gaussian")
                    plot.default (unlist(switched[[2]]), main="Residual Variance",
                                  xlab="Iteration",
                                  ylab="Residual Variance")
                main.label <- "Log-likelihood"; if (!is.null(DIC)) {
                                                    main.label <- paste0(main.label, ": DIC = ",signif(DIC[1], 5))
                                                    if (length(DIC)>=2) main.label <- paste0(main.label, "\nDeviance of Average = ",signif(DIC[2], 5))
                                                    if (length(DIC)>=3) main.label <- paste0(main.label, "\nEffective Parameter Count = ",signif(DIC[3], 5))
                                                }

                ## 3: Log-likelihood
                if (3 %in% which.plots) {
                    plot.default (unlist(switched[[3]]), main=main.label,
                                  xlab="Iteration",ylab="Log Likelihood",
                                  type="l")
                }

                ## 4: Ordinal cutoffs.
                if (4 %in% which.plots & class.outcome=="ordinal" & ordinal.count>2) {
                    draws <- length(unlist(switched[[1]]))
                    xx <- sort(rep(1:draws, ordinal.count-2))
                    plot.default (c(1,xx), c(0,unlist(switched[[4]])), col=c(0, rep(1:(ordinal.count-2), draws)),
                                  main="Ordinal Cutoff Values",
                                  xlab="Iteration",
                                  ylab="Cutoff Value")
                    abline(h=0, col=8)
                }

                ## 5: Reciprocity
                if (5 %in% which.plots & reciprocity.present) {
                    main.label <- "Reciprocal Correlation"
                    plot.default (unlist(switched[[5]]), main=main.label,
                                  xlab="Iteration",
                                  ylab="Reciprocal Correlation")


                }

                                        #6 to (5 + length(components)): components!
                if (length(components) > 0) for (cc in 1:length(components)) if ((non.comp.count()+cc) %in% which.plots) components[[cc]]$gibbs.plot(switched[[cc+non.comp.count()]])

            },


            DIC = function (gibbs.out, add.parts=FALSE) {
                                        #all.values <- gibbs.value(gibbs.out)
                                        #deviance.of.average <- -2*log.likelihood.by.value (apply(all.values, 1, mean))
                deviance.of.average <- -2 * gibbs.mean(gibbs.out)$log.likelihood
                average.deviance <- -2 * mean(unlist(gibbs.switcheroo(gibbs.out)$log.lik))
                                        #2*apply(all.values, 2, log.likelihood.by.value))
                output <- c(DIC=2*average.deviance - deviance.of.average)
                if (add.parts) output <- c(output,
                                           deviance.of.average=deviance.of.average,
                                           effective.parameters=average.deviance - deviance.of.average,
                                           average.deviance=average.deviance)
                return(output)
            },

            marginal.loglikelihood = function (gibbs.out) {
                all.values <- gibbs.value(gibbs.out)
                model.log.likelihoods <- apply(all.values, 2, log.likelihood.by.value)
                1/mean(1/model.log.likelihoods)
            },

            pseudo.CV.loglikelihood = function (gibbs.out) {
                all.values <- gibbs.value(gibbs.out)
                model.log.likelihoods <- apply(all.values, 2, log.likelihood.by.value, sumup=FALSE)
                each.like <- log(1/apply(1/exp(model.log.likelihoods), 1, mean))
                sum(each.like)
            }

        )
    )

### print.CIDnetwork <- function (x, ...) {}
### plot.CIDnetwork <- function (x, ...) {}
### summary.CIDnetwork <- function (object, ...) {}




###############################################################
####################  USER FUNCTIONS  #########################
###############################################################


                                        #CID <- function (...) CIDnetwork$new(...)
CID.generate <- function (...) CID (..., generate=TRUE)


CID <- function (input,    # Must be edge.list, sociomatrix, CID.object or CID.Gibbs.object
                 outcome,  # Only used when input is an edge.list
                 n.nodes,
                 node.names,
                 intercept = 0,# Only used if input is missing.
                 components,
                 class.outcome="ordinal",
                 fill.in.missing.edges=missing(outcome),
                 generate=FALSE,
                 verbose=2,
                 ...) {
                                        #  browser()
    if (missing(components)) components <- list()

    if (missing(input) & missing(n.nodes)) stop ("CID: no input was provided. Requires at least one of: sociomatrix, edge list, n.nodes.")
    if (!missing(input)){
        if(!(is(input,"matrix") | is(input,"data.frame") | is(input,"array"))) stop ("CID: input must be an edge.list or sociomatrix.")
    }
    if (missing(input)) {
        return(CIDnetwork$new(n.nodes=n.nodes, intercept=intercept,components=components, class.outcome=class.outcome, generate=TRUE, verbose=verbose, ...))
    } else { #get started!

        if ((is(input,"matrix") | is(input,"data.frame") | is(input,"array")) && ncol(input) == 2) { ##if we have an edgelist
            edge.list <- input

                                        #this is an impossible state      if (!(class(edge.list) %in% c("matrix", "array", "data.frame"))) stop ("The edge.list provided must be an array, matrix or data frame with two columns; this is not an array, matrix or data frame.")
                                        #also impossible      if (ncol(edge.list) != 2) stop ("The edge.list provided must have two columns.")

            edge.list <- cbind(as.character(edge.list[,1]), as.character(edge.list[,2]))
            if (any(edge.list[,1] == edge.list[,2])) {
                if(verbose > 0){
                    message ("Removing self edges.")
                }
                nonselfies <- which(edge.list[,1] != edge.list[,2])
                edge.list <- edge.list[nonselfies,]
                if (!missing(outcome)) outcome <- outcome[nonselfies]
            }

            node.names.tmp <- unique(c(as.character(edge.list[,1]), as.character(edge.list[,2])))
            is.num.list <- FALSE
            if(!missing(node.names)) {
                if(length(node.names) >= length(node.names.tmp)){
                    if(any(is.na(match(node.names.tmp,node.names)))){ ##node names don't match
                        suppressWarnings(node.nums <- as.numeric(node.names.tmp))
                        is.num.list <- !any(is.na(node.nums))
                        if(!is.num.list){  ## If edge.list has non-numbers, use edge.list names
                            if(verbose > 0)
                                message("Warning:  edge.list contains nodes not named in node.names.")
                            node.names <- node.names.tmp
                        }else{
                            node.names.tmp <- as.character(sort(node.nums))  ##  sorting the node name numbers
                            if(max(node.nums) > length(node.names)){
                                if(verbose > 0)
                                    message("Warning:  edge.list contains node numbers larger than node.names vector.")
                                node.names <- node.names.tmp
                            }
                        }
                    }
                }else{
                    if(verbose > 0)
                        message("Warning:  Length of node.names less than n.nodes.  Using default node names")
                    node.names <- node.names.tmp
                }
            }else{
                node.names <- node.names.tmp
            }

            n.nodes <- length(node.names)
            if(is.num.list){
                numbered.edge.list <- cbind(match(edge.list[,1], node.names.tmp),
                                            match(edge.list[,2], node.names.tmp))
            }else{
                numbered.edge.list <- cbind(match(edge.list[,1], node.names),
                                            match(edge.list[,2], node.names))
            }

            if (missing(outcome) | fill.in.missing.edges) {
                if(verbose > 0){
                    if (missing(outcome))
                        message("Assuming that this is a complete network with specified edges as binary ties.")
                    else
                        message("Filling in unspecified edges as zeroes.")
                }

                ##  OLD CODE THAT ASSUMED UNDIRECTED NETWORK
###        new.edge.list <- make.edge.list(n.nodes)
###        rowmatch <- sapply (1:nrow(edge.list),
###                            function(rr) min(which((numbered.edge.list[rr,1] == new.edge.list[,1] &
###                                                    numbered.edge.list[rr,2] == new.edge.list[,2]) |
###                                                   (numbered.edge.list[rr,2] == new.edge.list[,1] &
###                                                    numbered.edge.list[rr,1] == new.edge.list[,2]))))

                new.edge.list <- make.edge.list(n.nodes)
                new.edge.list <- rbind(new.edge.list,new.edge.list[,2:1])
                new.edge.list <- new.edge.list[order(new.edge.list[,1],new.edge.list[,2]),]


                rowmatch <- sapply (1:nrow(edge.list),
                                    function(rr) which((numbered.edge.list[rr,1] == new.edge.list[,1] &
                                                        numbered.edge.list[rr,2] == new.edge.list[,2])))

                rowmatch <- rowmatch[is.finite(rowmatch)]

                temp.outcome <- rep(0, nrow(new.edge.list))
                if (missing(outcome)) {
                    temp.outcome[rowmatch] <- 1
                } else {
                    temp.outcome[rowmatch] <- outcome[is.finite(rowmatch)]
                }

                outcome <- temp.outcome
                edge.list <- new.edge.list

            } else {
                edge.list <- numbered.edge.list
            }
        } else { ## we have a sociomatrix
            sociomatrix <- input
            if (nrow(sociomatrix) != ncol(sociomatrix)) stop ("The provided sociomatrix is not square.")

            n.nodes <- nrow(sociomatrix)

            if (any(sociomatrix[!is.na(sociomatrix)] != t(sociomatrix)[!is.na(t(sociomatrix))])) {
                if(verbose > 0)
                    message ("CID: Auto-detected an asymmetric sociomatrix; assuming this is a directed network.")
                edge.list <- make.arc.list(n.nodes)
                outcome <- t(sociomatrix)[non.diag(n.nodes)]
            } else {
                if(verbose > 0)
                    message ("CID: Auto-detected a symmetric sociomatrix; assuming this is an undirected network.")
                edge.list <- make.edge.list (n.nodes)
                outcome <- sociomatrix[u.diag(n.nodes)]  #just to be clear.
            }

            ## Inferring Node Names
            if(missing(node.names) || length(node.names) != ncol(sociomatrix)) {
                if(!missing(node.names)){
                    if(verbose > 0) message("Warning: Length of node.names differs from nodes in sociomatrix.  Using default node.names")
                }
                if(!is.null(colnames(sociomatrix))){
                    node.names <- colnames(sociomatrix)
                }else{
                    node.names <- 1:n.nodes
                }
            }
        }

        ordinal.count <- max(outcome,na.rm=TRUE)+1

        ##    if (is.null(class.outcome)) {

        if (class.outcome=="ordinal") {
            outcome.na <- outcome[!is.na(outcome)]
            if (!all(round(outcome.na)==outcome.na)) {
                if(verbose > 0)
                    message ("Detected non-integer values for outcomes; fitting this model with Gaussian outcomes.")
                class.outcome <- "gaussian"
            } else {
                if(verbose > 1)
                    message ("Fitting: ordinal outcome with ",ordinal.count," states.")
            }
        } else {
            if(verbose > 1)
                message ("Fitting: Gaussian outcome.")
        }

        CID.object <- CIDnetwork$new(edge.list, n.nodes=n.nodes,
                                     intercept=intercept,
                                     components=components, outcome=outcome,
                                     class.outcome=class.outcome,
                                     ordinal.count=ordinal.count,
                                        #                                 ordinal.cutoffs=ordinal.cutoffs,
                                     node.names=as.character(node.names),
                                     generate=generate,
                                     verbose=verbose, ...)

        return(CID.object)
    }

}


CID.Gibbs <- function (input, # Must be edge.list, sociomatrix, CID.object or
                                        # CID.Gibbs.object
                       outcome, # Only used when input is an edge.list
                       node.names,
                       ##edge.list,
                       ##outcome,
                       ##sociomatrix,

                       ##CID.object,
                       ##CID.Gibbs.object,

                       ##n.nodes=max(edge.list),

                       components, #=list(),
                       class.outcome="ordinal",
                       fill.in.missing.edges=missing(outcome),
                       new.chain=FALSE,

                       draws=100,
                       burnin=-1,  ## burnin = -1 performs auto-burnin.
                       thin=10,
                       report=100,
                       auto.converge=FALSE,
                       extend.max = 10,
                       extend.count = 100,
                       verbose = 2,

                       ...) {
                                        #edge.list=n0$edge.list; outcome=n0$outcome; components=list(); n.nodes=max(edge.list)
                                        #edge.list=dolphins; components=list(LSM(2)); class.outcome=NULL
                                        #data(prison); sociomatrix=prison

    if(!is(input,"CID.Gibbs")) {  #missing(CID.Gibbs.object)){
        if(!is(input,"CIDnetwork")) {  #then it's a sociomatrix/edge.list/outcome combo.
            CID.object <- CID (input,
                               outcome,
                               components=components,
                               class.outcome=class.outcome, #
                               node.names=node.names,
                               fill.in.missing.edges=fill.in.missing.edges,
                               verbose=verbose, ...)
        } else {
            CID.object <- input
            if (!missing(components)) CID.object$components <- components
            CID.object$reinitialize()
        }
        res <- CID.object$gibbs.full(draws=draws, burnin=burnin, thin=thin,
                                     report=report,
                                     auto.converge=auto.converge,
                                     extend.max = extend.max,
                                     extend.count = extend.count,
                                     verbose=verbose, ...)
    } else {
        if (!missing(components)) warning ("You cannot change the components of a CID.Gibbs object. Call this using the appropriate CID object if you wish to change the components.")
        CID.Gibbs.object <- input
        CID.object <- CID.Gibbs.object$CID.object
        res.2 <- CID.object$gibbs.full(draws=draws, burnin=burnin, thin=thin,
                                       report=report,
                                       auto.converge=auto.converge,
                                       extend.max = extend.max,
                                       extend.count = extend.count, ...)
        ##  Concatenating Chains
        if(new.chain){
            res <- res.2
        } else {
            res <- c(CID.Gibbs.object$results,res.2)
        }
    }

                                        # At this point we should have a CID object. Hurrah!


    DIC <- CID.object$DIC(res, add.parts=TRUE)

                                        #marg.ll <- CID.object$marginal.loglikelihood(res)
                                        #pseudo.CV.ll <- CID.object$pseudo.CV.loglikelihood(res)

    output <- list(results=res,  #unwrap.CID.Gibbs(
                   CID.object=CID.object,
                   CID.mean=CID.object$gibbs.mean(res),
                   DIC=DIC)
    class(output) <- "CID.Gibbs"

    return(output)
}


print.CID.Gibbs <- function (x, ...) {
                                        #x$CID.object$gibbs.summary(x$results, ...)
    x$CID.object$show(...)
}

summary.CID.Gibbs <- function (object, ...) {
    object$CID.object$gibbs.summary(object, ...)

}

print.summary.CID.Gibbs <- function(x, ...){
    x$CID.object$print.gibbs.summary(x)
}


plot.CID.Gibbs <- function (x, ...) {
    x$CID.object$gibbs.plot (x$results, x$DIC, ...)
}


## predict.CID.Gibbs <- function(x, ...) {

## }