.packageName <- "gRbase"
  ## ####################################################################
  ## Interface   gRbase <-> dynamicGraph
  ## ####################################################################

UserMenus <- list(MainUser =
                  list(label = "Fit",
                       command = function(object, ...)
                       {
                         args <- list(...)
                         env <- args$Arguments
                         
                         object.new <- fit(object)
                         
                         env$redrawGraphWindow(
                                               env$graphLattice,
                                               env$graphWindow,
                                               env$edgeList,
                                               env$blockEdgeList,
                                               env$factorVertexList,
                                               env$factorEdgeList,
                                               env$visibleVertices,
                                               env$extraList,
                                               object = object.new,
                                               title = "Result from Fit",
                                               objectName = "currentObject",
                                               Transformation = NULL,  
                                               background = "white",
                                               vertexcolor = "black",
                                               w = 10, width = 400,  
                                               height = 400)                                     
                         return(list(object=object.new))
                       },
                       update.vertices = TRUE,
                       update.edges = TRUE
                       ),
                  MainUser = list(label = "Label all edges", 
           command = function(object, ...) LabelAllEdges(object, 
                    slave = FALSE, ...)),

                  )



gRVariableDescription <- function(obj) {
  ## obj is an hllm object
  
  object <- gmData(obj)
  
  return(list(
              names = as.character(varNames(object)),
              labels = NULL,
              types = as.character(varTypes(object)),
              levels = valueLabels(object)
              )
         )
}

gREdges <- function(object)
  {
    ## object is an hllmclass
    nodelabels <- varNames(gmData(object))
    form <- Formula(object)
    listform <- readf(form[2])
#              new.form <- add.edge(listform,c(name.1,name.2))
#              print(listform)
    from <- c()
    to <- c()
    for (i in 1:length(listform)) {
#                print(listform[[i]])
      edges <- selectOrder(listform[[i]])
      edges.ul <- unlist(edges)
#      print(edges.ul) 
      from <- c(from,edges.ul[1:length(edges.ul)%%2==1])
#      print(1:length(edges.ul)%%2)
#      print((1:length(edges.ul)+1)%%2)
      to <- c(to,edges.ul[(1:length(edges.ul)+1)%%2==1])
#      print(c(1:length(nodelabels))[unlist(edges)==nodelabels])
    }
#    object.UG <- as.UG(Formula(object))
#    edgemat <- allEdges(object.UG)


#      print(from)
#      print(to)
    
    from <- match(from,nodelabels)
    to   <- match(to,nodelabels)
    edgemat <- cbind(from,to)

#      print(from)
#      print(to)

#    print(edgemat)
    
#    labels.gR <- varNames(gmData(object))    
#    labels <- rownames(object.UG)

    
    ## convert ggm-edges to gR-edges
#    ggm2gr <- function(edge) match(labels[edge],labels.gR)
    
    if (length(edgemat)==0) return(matrix(nrow=0,ncol=2))
    
#    return(t(apply(edgemat,1,ggm2gr)))
#    print(edgemat)
        return(edgemat)
  }

dynamic.gR.Graph <-  function(object, ...)
{

  require(dynamicGraph)
  .Load.gRbase.dynamic()
    
    if (inherits(object,"gmData")) 
      object <- new("hllm",~.^.,object)

  VariableDescription <- gRVariableDescription(obj = object)

   if (!is.graphical(readf(Formula(object)[2]))) {
     cat("Model not graphical, using factorgraph\n")
     ## factor-edges, factorgraph
#     FactorEdges <- gRFactorEdges(object = object)
     Edges<-NULL
     FactorEdges <-     readf(Formula(object)[2])
   }
  else {
    Edges <- gREdges(object = object)
    FactorEdges <- NULL
  }

#    print(VariableDescription)
#    print(Edges)
#    print(FactorEdges)
    
    Z <- DynamicGraph(names = VariableDescription$names,
                      types = VariableDescription$types,
                      from = Edges[,1], to = Edges[,2],
                      factors=FactorEdges,
                      oriented = FALSE, 
                      object = object,
                      UserMenus = UserMenus,
                      ...)
  }

LabelAllEdges <- function(object, slave = FALSE, 
                ...) {
                args <- list(...)
                Args <- args$Arguments
                getNodeName <- function(index, type) if (type == 
                  "Vertex") 
                  name(Args$vertexList[[index]])
                else if (type == "Factor") 
                  name(Args$factorVertexList[[abs(index)]])
                else if (type == "Block") 
                  label(Args$blockList[[abs(index)]])
                else NULL
                visitEdges <- function(edges) {
                  for (i in seq(along = edges)) {
                    vertices <- nodeIndicesOfEdge(edges[[i]])
                    types <- nodeTypesOfEdge(edges[[i]])
                    name.f <- getNodeName(vertices[1], types[1])
                    name.t <- getNodeName(vertices[2], types[2])
                    R <- testEdge(object, action = "remove", 
                      name.1 = name.f, name.2 = name.t, from = vertices[1], 
                      to = vertices[2], from.type = types[1], 
                      to.type = types[2], edge.index = i, force = force, 
                      Arguments = Args)
                    if (!is.null(R)) {
                      if (TRUE || (hasMethod("label", class(R)))) 
                        label(edges[[i]]) <- label(R)
                      if (TRUE || (hasMethod("width", class(R)))) 
                        width(edges[[i]]) <- width(R)
                    }
                  }
                  return(edges)
                }
                edgeList <- visitEdges(Args$edgeList)
                factorEdgeList <- visitEdges(Args$factorEdgeList)
                blockEdgeList <- visitEdges(Args$blockEdgeList)
                if (slave) 
                  Args$redrawGraphWindow(graphWindow = NULL, 
                    edgeList = edgeList, factorEdgeList = factorEdgeList, 
                    blockEdgeList = blockEdgeList, title = "A slave window", 
                    ...)
                else Args$redrawGraphWindow(graphWindow = Args$graphWindow, 
                  edgeList = edgeList, factorEdgeList = factorEdgeList, 
                  blockEdgeList = blockEdgeList, title = "Not used!", 
                  width = NULL, height = NULL, Arguments = Args)
            }
                
## ##########################################################
## ###### formulae #############
## ##########################################################
  
## Interface gRbase <-> ggm
##
## This should be replaced with an interface to gRaph, since ggm
## uses adjancancy-matrix representation and not lists as gRaph does.
## ####################################################################


all.subsets <- function(x,g.sep="+"){
  if (length(x)==1)
    return(x)
  else {
    val <- x[1]
    for (i in 2:length(x)){
      v <- paste(val,x[i],sep=g.sep)
      val <- c(val,x[i],v)
    }
    val <- strsplit(val,paste("\\",g.sep,sep=""))
    return(val)
  }
}

selectOrder  <- function(x,order=2){
  v <- all.subsets(x)
  ##print(x); print(v); print(order)
  value <- v[lapply(v,length)==as.numeric(order)]
  return(value)
}

extract.power<-function(fff){
  mimf  <- paste(as.formula(fff))[2]
  mimf.split <- unlist(strsplit(mimf,""))
  if(length(grep("[a-z]", mimf))>0){
    pow <- mimf 
  } else {
    has.hat <- match("^",mimf.split)
    sub <- unlist(strsplit(mimf,"\\^"))
    ##print(sub)
    if (!is.na(has.hat)){
      pow <- ifelse (sub[2]==".", -1, as.numeric(sub[2]))
    }
    else {
      pow <- length(unlist(strsplit(sub,"\\.")))
    }
  }
  return(pow)
}



process.formula <- function(formula, data, marginal, type=c("Discrete","Continuous"),v.sep=":",g.sep="+"){
  
  get.var.of.type <- function(type){varNames(data)[varTypes(data)==type]}
  
  used.var <- get.var.of.type(type)
  pow <- extract.power(formula)
  if (is.numeric(pow)){
    if (!missing(marginal)){
      ##      print(used.var); print(marginal)
      used.var <- intersect(marginal,used.var)
    }
    if (pow==-1)
      mimf <- paste(used.var,collapse=v.sep,sep="")
    else{
      pow <- min(c(pow, length(used.var)))
      tmp <- selectOrder(used.var, pow)
      mimf <- paste(unlist(lapply(tmp, paste, collapse=v.sep)),collapse=g.sep,sep="")
    }
  } else {
    mf    <- as.formula(formula)
    mimf  <- paste(mf,sep="")[2]
  }
  formula <- formula(paste("~",mimf,sep=""))
  interactions <- strsplit(mimf,paste("\\",g.sep,sep=""))[[1]]
  interactions <- gsub(g.sep,"",interactions)
  int.list <- strsplit(interactions, v.sep)
  gc1   <- lapply(int.list, function(l){ match(l,used.var) })
  
  value <- list(formula=formula, mimformula=mimf, numformula=gc1,
                gmData=data, varnames=used.var)
  value
}

## ##########################################################
## ###### end formulae #############
## ##########################################################

# Notes, functions and examples for generator lists (model formulae
# for hierarchical loglinear models) in R. 
# Includes functions dual.rep, add.edge, delete.edge and is.graphical.
# I havent tried to optimise the functions in any way, merely to get
# versions which work. 
# David Edwards, 12.5.2004.
# adapted to gRbase by Claus Dethlefsen 16.08.04

# The approach could also be used (with some extra work) for
# hierarchical mixed 'mim' models,  
# but with 'extended' hierarchical mixed models there is a problem of
# how to handle quadratic terms like X^2. 
# Representing such generators as vectors would seem problematic.

# nb: there is an implicit assumption that all models contain at least
# all main effects, eg. the minimal model 
# for variable set A,B,C has formula A+B+C.

# ----------------------------------------------------------
# ---------------------------------------------------------- 
#                         Generators
# ----------------------------------------------------------
# ----------------------------------------------------------          
 
# A generator is implemented as a vector (regarded as a set)
#
#g1 <- c("a", "bc", "x")
#g2 <- c("a", "x")
#g3 <- c("sex", "Age")
#g4 <- 1:5
#
# Thus the following built-in R commands may be used with generators:
#
#    union(g1, g2)
#    intersect(g1, g2)
#    setdiff(g1, g2)
#    setequal(g1, g2)
#    is.element(g1, g2)
#
# setdiff(g1,g2) returns g1 \ g2.
# Comparing two generators, is.element(g1, g2) returns a boolean 
# vector of the same length as g1, indicating whether the element of
# g1 is contained in g2. Thus
# g1 <= g2 <=> all(is.element(g1, g2))
# g1 == g2 <=> setequal(g1, g2)
#
subsetof <- function(g1, g2) all(is.element(g1, g2))  
  
# A function to write a generator as a string:
#
showg <- function(g, v.sep=":") {
  if (length(g)==0)
    s<-'<empty>'
  else {s <- g[1];
        if (length(g)>1)
          for (i in 2:length(g))
            s <- paste(s, g[i], sep=v.sep)
      }
  s
}

# Sometimes we may need the empty set
#g5 <- vector()
#showg(g5)

# A function to read a generator as a string. 
# nb: s is a character (vector), but must have length one.
#
readg <- function(s, v.sep=":") {
  g <- character(0)
  s <- paste(s, v.sep, sep="") # add a separator
  s <- gsub(" ","",s) # strip spaces
  k1 <- 1
  for (k in ((k1+1):nchar(s))) {
    if (substring(s,k,k) == v.sep) {
      g <- c(g, substring(s,k1,(k-1)))
      k1 <- k+1 }    
  }
  g
}

# -------------------------------------------------------------------------------------------------------
#               Generator Lists
# -------------------------------------------------------------------------------------------------------

#f1 <- list(g1, g2, g3)
#f2 <- list(1:10, c(2,3,5), c(3,5,7))


# Function to write a generator list as a formula
showf <- function(f, g.sep="+", v.sep=":") {
   if (length(f)==0) s <- '<empty formula>' else {
   s <- showg(f[[1]])
   if (length(f)>1)
     for (i in 2:length(f))
       s <- paste(s, showg(f[[i]],v.sep), sep=g.sep)}
   s
}

# Function to read a generator list as a string 
# nb: length of string must be one

readf <- function(s, v.sep=":", g.sep="+") {
  gens <- readg(s, g.sep) 
  l <- list(length(gens))
  for (i in 1:length(gens)) l[[i]] <- readg(gens[i], v.sep)
  l
}

#f3 <- readf("A.B.C+B.C.D+C.D.E")
#showf(f3)

# To get the union of all generators i a list l:

varset <- function(f) unique(unlist(f))

# To find out whether a generator g is contained in (is a subset of
# some element of) a list l: 

in.list <- function(g, l)
  any(unlist(lapply(l, function(xx) all(is.element(g, xx)))))

# To find out whether a generator g contains an element of a list l:

# any(unlist(lapply(l, function(xx) any(is.element(xx, g))))

# A function to find out whether the k.th generator in a list 
# is contained in any of the others: 

is.cont <- function(k, l) {
  g <- l[[k]]
  a <- sapply(l, function(xx) all(is.element(g, xx)))
  a[k] <- F
  any(a)
}

# A function to find out whether the k.th generator in a list
# contains any of the others 

contains <- function(k, l) {
  g <- l[[k]]
  a <- unlist(lapply(l, function(xx) all(is.element(xx, g))))
  a[k] <- F
  any(a)
}

# A function to remove redundant generators.
# If maximal=T, returns the maximal generators, if =F, the minimal generators.
# This seems difficult to vectorize: any suggestions? 
# This function is a prime candidate for a C routine.

remove.redundant <- function(f, maximal=TRUE) {
  k <- length(f)
  new.f <- f
  if (k>1) {
    for (i in 1:(k-1)) {
      g1 <- f[[i]]
      if (length(g1)>0) {
        for (j in (i+1):k) {
          g2 <- f[[j]]
          if (length(g2)>0) {
           if (setequal(g1,g2)) f[[j]]<-vector() else {
             kk <- 0
             if (subsetof(g1, g2)) {if (maximal) kk <- i else kk <- j}
             if (subsetof(g2, g1)) {if (maximal) kk <- j else kk <- i}
             if (kk>0) f[[kk]] <- vector()
      } # else    
     }} # length(g2)>0; for j in (i+1):k
    }}  # length(g1)>0; for i in 1:(k-1)
    f <- f[lapply(f, length)>0]
  } # if k>1
  f
}

# A function to return dual representation for a list. 
# See description in Edwards & Havranek, Biometrika (1985), 72, 2, p.341.

# minimal=T: returns dual representation given usual, ie list of
# minimal generators not contained in a generator in the input list. 
# minimal=F: returns usual representation given dual, ie list of
# maximal generators not containing a generator in the input list. 

## old version gave list()
#dual.rep <- function(glist, S, minimal=TRUE) {
# # S is total varset - often but by no means always given by
# # unique(unlist(g.list))  
# list.save <- list()
# if (length(glist)>0) for (v in 1:length(glist)) {
#   m1 <- list.save
#   if (minimal) m2 <- as.list(setdiff(S,glist[[v]])) else m2 <- as.list(glist[[v]])
#   if (v==1) list.save <- m2 else {
#      list.save <- remove.redundant(unlist(lapply(m1, function(g) lapply(m2, union, g)),recursive=FALSE),FALSE)}}
# if (!minimal) list.save <- lapply(list.save, function(g) setdiff(S, g))
# list.save
#} 

dual.rep <- function(glist, S, minimal=TRUE) {
 # S is total varset - often but by no means always given by unique(unlist(g.list)) 
 list.save <- list()
 #if (length(glist)==0) list.save <- list(S)
 if (length(glist)==1 & is.logical(glist[[1]])) list.save <- list(S)
 else { 
   for (v in 1:length(glist)) {
     m1 <- list.save
   if (minimal) m2 <- as.list(setdiff(S,glist[[v]])) else m2 <- as.list(glist[[v]])
   if (v==1) list.save <- m2 else {
      list.save <- remove.redundant(unlist(lapply(m1, function(g)
                                                  lapply(m2, union,
                                                         g)),recursive=FALSE),FALSE)}}
 if (!minimal) list.save <- lapply(list.save, function(g) setdiff(S,
                                                                  g))}
 list.save }  


#m1 <- readf('A.B+A.C')
#showf(m1)
#showf(dual.rep(m1, varset(m1)))

#m2 <- readf('B.D+A.D+C.D')
#showf(m2)
#showf(dual.rep(m2, varset(m2)))

# The dual of the dual should be the same as the original
#showf(dual.rep(dual.rep(m2, varset(m2)), varset(m2), F))

# Function to delete 'edge' from a generator list, by (i) converting
# generator list to dual representation, 
# (ii) appending the 'edge', (iii) converting back to usual
# representation. 'Edge' is given as vector of length 2:  
# it can also have length >2, ie. be a higher-order interaction.

delete.edge <- function(m, edge) {
  S <- varset(m)
  dr <- dual.rep(m, S)
  dr <- c(dr, list(edge))
  remove.redundant(dual.rep(dr, S, FALSE)) 
}

#m2 <- readf('B.D+A.D+C.D')
#showf(m2)
#showf(delete.edge(m2, c('A','B')))
#showf(delete.edge(m2, c('A','D')))
#m3 <- readf('A.B.C.D')
#showf(delete.edge(m3, c('A','B','C')))

# Function to add 'edge' from a generator list, by (i) converting
# generator list to dual representation, 
# (ii) removing the 'edge', (iii) converting back to usual
# representation. 'Edge' is given as vector of length 2:  
# it can also have length >2, ie. be a higher-order interaction.

add.edge <- function(m, edge) {
  S <- varset(m)
  dr <- dual.rep(m, S)
  k <- length(dr)
  if (k>0) {for (i in 1:k) if (setequal(dr[[i]], edge)) dr[[i]] <- vector()} 
#  if (k>0) {for (i in 1:k) if (setequal(dr[[i]], edge)) dr[[i]] <- NULL}
  dr <- remove.redundant(dr, FALSE)
  dual.rep(dr, S, FALSE)
}

#m2 <- readf('B.D+A.D+C.D')
#showf(m2)
#showf(add.edge(m2, c('A','B')))
#showf(add.edge(m2, c('A','C')))


# From main effect model on 5 vertices
#showf(add.edge(1:5, c(1,2)))
#showf(add.edge(c('a','b','c','d','e'), c('a','b')))

# From saturated model on 5 vertices
#showf(delete.edge(list(1:5), c(1,2)))
#showf(add.edge(list(c('a','b','c','d','e')), c('a','b')))

# Exploiting the fact that a model is graphical iff all its dual
# generators have length 2, 
# we get a neat function for graphicalness:
 
is.graphical <- function(m) {
   dr <- dual.rep(m, varset(m))
   lengths <- lapply(dr, length)
   all(lengths==2)
}

#is.graphical(readf('A.B+A.C+B.C'))
#is.graphical(readf('A.B.C'))
#is.graphical(readf('A.B+B.C'))
.Load.gRbase.general <- function() {

  ## ##########################################################
  ## ###### gmData #############
  ## ##########################################################
  
  ## a virtual class for data in any format (or NULL)
  setClassUnion("dataOrNULL", c("NULL","data.frame","table"))
  
  ##setClassUnion("displayOrNULL", c("NULL","dynamicDisplay"))
  
  ## Other datatypes may be added later:
  ## setIs("mydataclass", "dataOrNULL")
  ##
  ## existing may be inspected as
  ## getClass("dataOrNULL")

  ## ####################################################################
  ## CLASS: gmData. 
  setClass('gmData', representation(description='data.frame',
                                    valueLabels='vector',
                                    observations='dataOrNULL'
                                    )
           )
  
  ## creator. Default: no data 
  setMethod("initialize", "gmData", function(.Object,
                                             varNames=vector(),
                                             varTypes=
                                             rep(validVarTypes()[1],
                                                 length(varNames)),
                                             numberLevels=NA,
                                             latent=FALSE,
                                             valueLabels=list(),
                                             observations=NULL) {
    
    .Object@description              <- data.frame(I(varNames))
    .Object@description$numberLevels <- numberLevels
    .Object@description$latent       <- latent
    .Object@description$varTypes <- factor(varTypes,levels=validVarTypes())
    
    ## TODO: create numberLevels=2 for discrete variables
    ## TODO: create valueLabels for discrete variables
    
    .Object@valueLabels   <- valueLabels
    .Object@observations  <- observations
    
    .Object
  }
            )

  ## note that 'description' can be extended by further information
  ## about the variables.
  ##
  ## description(mydata)$mimName <- letters[1:nrow(description(mydata))]
  
  ## ####################################################################
  ## get/set methods
  
  ## description

  if (!isGeneric("description")) {
    if (is.function("description")) 
      fun <- description
    else fun <- function(x) standardGeneric("description")
    setGeneric("description", fun)
  }
  setMethod("description","gmData",function(x) x@description)
  
  setGeneric("description<-", function(x, value) standardGeneric("description<-"))
  setReplaceMethod("description", "gmData", function(x, value){
    x@description <- value
    x
  })
  
  ## get/set variable properties

  ## varTypes
  
  if (!isGeneric("varTypes")) {
    if (is.function("varTypes")) 
      fun <- varTypes
    else fun <- function(x) standardGeneric("varTypes")
    setGeneric("varTypes", fun)
  }
  setMethod("varTypes","gmData",function(x) description(x)$varTypes)
  
  setGeneric("varTypes<-", function(x, value) standardGeneric("varTypes<-"))
  setReplaceMethod("varTypes", "gmData", function(x, value){
    tmp <- 
      description(x)$varTypes <- value
    x
  })
  
  ## varNames
  
  if (!isGeneric("varNames")) {
    if (is.function("varNames")) 
      fun <- varNames
    else fun <- function(x) standardGeneric("varNames")
    setGeneric("varNames", fun)
  }
  setMethod("varNames","gmData",function(x) description(x)$varNames)
  
  setGeneric("varNames<-", function(x, value) standardGeneric("varNames<-"))
  setReplaceMethod("varNames", "gmData", function(x, value){
    description(x)$varNames <- value
    x
  })
  
  ## numberLevels
  
  if (!isGeneric("numberLevels")) {
    if (is.function("numberLevels")) 
      fun <- numberLevels
    else fun <- function(x) standardGeneric("numberLevels")
    setGeneric("numberLevels", fun)
  }
  setMethod("numberLevels","gmData",function(x) description(x)$numberLevels)
  
  setGeneric("numberLevels<-", function(x, value) standardGeneric("numberLevels<-"))
  setReplaceMethod("numberLevels", "gmData", function(x, value){
    description(x)$numberLevels <- value
    x
  })
  
  ## latent
  
  if (!isGeneric("latent")) {
    if (is.function("latent")) 
      fun <- latent
    else fun <- function(x) standardGeneric("latent")
    setGeneric("latent", fun)
  }
  setMethod("latent","gmData",function(x) description(x)$latent)
  
  setGeneric("latent<-", function(x, value) standardGeneric("latent<-"))
  setReplaceMethod("latent", "gmData", function(x, value){
    description(x)$latent <- value
    x
  })
  
  ## valueLabels
  
  if (!isGeneric("valueLabels")) {
    if (is.function("valueLabels")) 
      fun <- valueLabels
    else fun <- function(x) standardGeneric("valueLabels")
    setGeneric("valueLabels", fun)
  }
  setMethod("valueLabels","gmData",function(x) x@valueLabels)
  
  setGeneric("valueLabels<-", function(x, value) standardGeneric("valueLabels<-"))
  setReplaceMethod("valueLabels", "gmData", function(x, value){
    x@valueLabels <- value
    x
  })

  
  ## observations
  
  if (!isGeneric("observations")) {
    if (is.function("observations")) 
      fun <- observations
    else fun <- function(x) standardGeneric("observations")
    setGeneric("observations", fun)
  }
  setMethod("observations","gmData",function(x) x@observations)
  
  setGeneric("observations<-", function(x, value) standardGeneric("observations<-"))
  setReplaceMethod("observations", "gmData", function(x, value){
    x@observations <- value
    
    ## How do we ensure that varNames that are not latent can be
    ## found in data (should we?)?
    
    x
  })
  
  
  
  setMethod("show","gmData", function(object) {
    cat("Description:\n")
    show(description(object))
    cat("To see the values of the factors use the 'valueLabels' function\n")
    cat("To see the data use the 'observations' function\n")
}
            )
  
  ## ####################################################################
  ## Convert data.frame into gmData
  
  setAs("data.frame","gmData", function(from,to) {
    
    fact   <- unlist(lapply(1:ncol(from), function(j)
                            is.factor(from[,j])))
    Types <- rep(validVarTypes()[3],length(fact))
    Types[fact] <- validVarTypes()[1]
    
    levels <- unlist(lapply(1:ncol(from),
                            function(j)
                            {
                              if(is.factor(from[,j]))
                                length(levels(from[,j]))
                              else NA}
                            )
                     )
    
    if (length(which(fact))>0){
      vallabels <- list()
      for (j in which(fact)){
        vallabels <- c(vallabels, list(levels(from[,j])))
      }
      names(vallabels) <- names(from[which(fact)])
    } else {
      vallabels <- list()
    }
    
    new("gmData",
        varNames=names(from),
        varTypes=Types,
        numberLevels=levels,
        valueLabels=vallabels,
        observations=from
        )
  }
        )
  
  ## ####################################################################
  ## Convert table into gmData
  
  setAs("table","gmData", function(from,to) {
    counts <- as.vector(from)
    dn     <- dimnames(from)
    name   <- names(lapply(dn,function(x)names(x)))
    dim    <- unlist(lapply(dn,length))
    new("gmData",
        varNames=name,
        varTypes="Discrete",
        numberLevels=dim,
        valueLabels=dn,
        observations=from
        )
  }
        )
  
  ## ####################################################################
  ## Convert array into gmData
  ## example of adding a data-type
  
  setIs("array", "dataOrNULL")
  
  setAs("array","gmData", function(from,to) {
    res <- as(as.table(from),"gmData")
    observations(res) <- from
    res
  }
        )
  
  ## ##########################################################
  ##  ###### end gmData #############
  ## ##########################################################
  
  ## ####################################################################
  ## S4 class for gModel -- graphical models
  ## ####################################################################
  setClass("gModel", representation(
                                    formula = "formula",
                                    gmData = "gmData"
                                    )
           )

  
  
  ## ####################################################################
  ## end S4 class for gModel -- graphical models
  ## ####################################################################
  
    ## get function - gmData
  if(!isGeneric("gmData")){
    if (is.function("gmData"))
      fun <- gmData
    else fun <- function(object) standardGeneric("gmData")
    setGeneric("gmData", fun)
  }
  

  ## get function - formula
  if(!isGeneric("Formula")){
    if (is.function("Formula"))
      fun <- Formula
    else fun <- function(object) standardGeneric("Formula")
    setGeneric("Formula", fun)
  }
  
  
  ## set function - formula
  
  setGeneric("Formula<-", function(x, value) standardGeneric("Formula<-"))



  
  ## #############################
  ## Modify methods 
  ## #############################
  
  if(!isGeneric("dropEdge")){
    if (is.function("dropEdge"))
      fun <- dropEdge
    else fun <- function(object,name.1,name.2) standardGeneric("dropEdge")
    setGeneric("dropEdge", fun)
  }
  
  
  
  if(!isGeneric("addEdge")){
    if (is.function("addEdge"))
      fun <- addEdge
    else fun <- function(object,name.1,name.2) standardGeneric("addEdge")
    setGeneric("addEdge", fun)
  }
  
  
  if(!isGeneric("dropVertex")){
    if (is.function("dropVertex"))
      fun <- dropVertex
    else fun <- function(object,name) standardGeneric("dropVertex")
    setGeneric("dropVertex", fun)
  }
  
  
  if(!isGeneric("addVertex")){
    if (is.function("addVertex"))
      fun <- addVertex
    else fun <- function(object,name) standardGeneric("addVertex")
    setGeneric("addVertex", fun)
  }
  
}

.Load.gRbase.dynamic <- function() {
  require(dynamicGraph)
  
  if (!isGeneric("label") && !isGeneric("label", where = 3)) {
    if (is.function("label"))
      fun <- label
    else
      fun <- function(object) standardGeneric("label")
    setGeneric("label", fun)
  }
  
  setMethod("label", "hllmTestClass", function(object)
            format(object@p, digits = 4))
  
  if (!isGeneric("width") && !isGeneric("width", where = 3)) {
    if (is.function("width"))
      fun <- width
    else
      fun <- function(object) standardGeneric("width")
    setGeneric("width", fun)
  }
  
  setMethod("width", "hllmTestClass", function(object)
            round(2 + 5 * (1 - object@p)))


  
}

.Load.dynamicgraph <- function() {
    if (!isGeneric("dynamic.Graph")) {
    if (is.function("dynamic.Graph")) 
      fun <- dynamic.Graph
    else fun <- function(object, ...) standardGeneric("dynamic.Graph")
    setGeneric("dynamic.Graph", fun)
  }
  setMethod("dynamic.Graph", signature(object = "hllm"), 
            function(object, ...) {
              dynamic.gR.Graph(object, title="Hierarchical log-linear model",...)
        })
  setMethod("dynamic.Graph", signature(object = "gmData"), 
            function(object, ...) {
              dynamic.gR.Graph(new("hllm",.^.~1,object), title="Hierarchical log-linear model",...)
        })

  }
.Load.gRbase.hllm <- function() {

    ## ####################################################################
  ## S4 class for hllm
  ## ####################################################################
  
  setClass("hllm", contains="gModel")
           
  ## creator
  setMethod("initialize", "hllm", function(.Object,formula=~.^1,gmData,marginal)
            {
              .Object@formula <- process.formula(formula,gmData,marginal,type="Discrete")$formula
              .Object@gmData  <- gmData
              .Object
            }
            )

    setMethod("gmData", "hllm", function(object) object@gmData)
    setMethod("Formula", "hllm", function(object) object@formula)
  setReplaceMethod("Formula", "hllm", function(x,value) {
    
    x@formula <- value
    x
  })

setMethod("dropEdge", "hllm",
          function(object,name.1,name.2) {
            
            ## cat("Drop:",name.1,name.2,"\n",sep=" ")
              
            ## edit hllm formula
            form <- Formula(object)
            listform <- readf(form[2])
            new.form <- delete.edge(listform,c(name.1,name.2))
              
            form <- paste("~",showf(new.form))
            Formula(object) <- as.formula(form)
            new.object <- object
              
            if (extends(class(new.object),"hllmengine"))
              new.object <- fit(new.object)
            
            return(new.object)
          }
          )

setMethod("addEdge", "hllm",
          function(object,name.1,name.2) {
            
            ## edit hllm formula
            form <- Formula(object)
            listform <- readf(form[2])
            new.form <- add.edge(listform,c(name.1,name.2))
            form <- paste("~",showf(new.form))
            Formula(object) <- as.formula(form)
            new.object <- object
            
            if (extends(class(new.object),"hllmengine"))
              new.object <- fit(new.object)
            return(new.object)
          })


setMethod("dropVertex", "hllm",
          function(object,name) {
            ## edit hllm formula
            form <- Formula(object)
            listform <- readf(form[2])

            ## delete 'name' from generators
            new.form <- lapply(listform,setdiff,name)
            form <- paste("~",showf(new.form))
            Formula(object) <- as.formula(form)
            
            new.object <- object
              
            if (extends(class(new.object),"hllmengine"))
              new.object <- fit(new.object)
            return(new.object)
          })

setMethod("addVertex", "hllm",
          function(object,name) {
            ## edit hllm formula
            form <- Formula(object)
            listform <- readf(form[2])
            listform[[length(listform)+1]] <- name
            form <- paste("~",showf(listform))
            Formula(object) <- as.formula(form)
#              u <- as.UG(form)
#              u <- cbind(u,0)
#              u <- rbind(u,0)
#              rownames(u)[nrow(u)] <- name
#              colnames(u)[ncol(u)] <- name
#              u.cl <- cliques(u)
#              u.form <- paste(unlist(lapply(u.cl,paste,collapse=":")),collapse=" + ")
#              form <- paste(u.form,form[1],form[3])
#            form <- paste("~",u.form)
#              Formula(object) <- as.formula(form)
#            object@formula <- as.formula(form)
              new.object <- object
              
              if (extends(class(new.object),"hllmengine"))
                new.object <- fit(new.object)
              return(new.object)
            })

  
  ## ####################################################################
  ## hllmTest  (adopted from CoCoObjects)
  ## ####################################################################

  setClass("hllmTestClass", representation(deviance = "numeric", 
                                           df = "numeric", p = "numeric"))

  ## ####################################################################
  ## end hllmTest  (adopted from CoCoObjects)
  ## ####################################################################

}


.Load.gRbase.hllmfit <- function() {
    ## ####################################################################
  ## hllm fit
  ## ####################################################################
  
  ## a virtual class for fit output in any format depending on the
  ## engine
  setClassUnion("hllmengine")
  
  if(!isGeneric("getFit")){
    if (is.function("getFit"))
      fun <- getFit
    else fun <- function(object,...) standardGeneric("getFit")
    setGeneric("getFit", fun)
  }
  
  setMethod("getFit","hllm",function(object) object@fit)

  if(!isGeneric("summary")){
    if (is.function("summary"))
      fun <- summary
    else fun <- function(object,...) standardGeneric("summary")
    setGeneric("summary", fun)
  }

setMethod("summary","hllm",function(object) summary(getFit(object)))

  
  ## test
  ## extends(class(currentObject),"hllmengine")
  
  
  if(!isGeneric("fit")){
    if (is.function("fit"))
      fun <- fit
    else fun <- function(object,...) standardGeneric("fit")
    setGeneric("fit", fun)
  }
  
  setMethod("fit", "hllm",
            function(object,engine="loglm",...) {
              obj <- as(object,paste("hllm",engine,sep=""))
              res <- fit(obj,...)
              res
            }
            )

  ## loglin engine
  setClass("hllmloglin", contains="hllm",representation(fit="list"))
  setIs("hllmloglin","hllmengine")
  setIs("hllmloglin","hllm")
  
  setAs("hllm","hllmloglin", function(from,to) {
    new("hllmloglin",Formula(from),gmData(from))
  }
        )
  setMethod("fit","hllmloglin", function(object,...) {
    
    rawdata <- observations(gmData(object))
    if (is.data.frame(rawdata)){
      rawdata <- xtabs(~., rawdata)
      ##rawdata <- as.data.frame(rawdata)
    }
    numform <- process.formula(Formula(object),gmData(object),type="Discrete")$numformula
    val <- loglin(rawdata, numform,...)
#    print(class(val))
    object@fit <- val
    object
  }
            )
  
  ## loglm engine
  setClass("hllmloglm", contains="hllm",representation(fit="loglm"))
  setIs("hllmloglm","hllmengine")
  setIs("hllmloglm","hllm")
  
  setAs("hllm","hllmloglm", function(from,to) {
    new("hllmloglm",Formula(from),gmData(from))
  }
        )
  setMethod("fit","hllmloglm", function(object,...) {
    require(MASS)
    
    rawdata <- observations(gmData(object))
    if (is.data.frame(rawdata)){
      rawdata <- xtabs(~., rawdata)
      ##rawdata <- as.data.frame(rawdata)
    }
    mimform <- process.formula(Formula(object),gmData(object),type="Discrete")$mimformula
    loglm.formula <- formula(paste("~",mimform))
    val <- loglm(loglm.formula, rawdata,...) 
    
    object@fit <- val
    object
  }
            )
  
  ## ####################################################################
  ## end hllm fit
  ## ####################################################################

}

.Load.gRbase.hllmodify <- function() {
    ## ####################################################################
  ## Editing hllm object 
  ## ####################################################################
  
  ## modifyModel is identical to the one from CoCo, except the
  ## class. Except also that subModifyModel is replaced with a method
  ## for each of the possibilities.
  if (!isGeneric("modifyModel")) {
    if (is.function("modifyModel")) 
      fun <- modifyModel
    else fun <- function(object, action, name, name.1, name.2, 
                         ...) standardGeneric("modifyModel")
    setGeneric("modifyModel", fun)
  }
  setMethod("modifyModel", signature(object = "hllm"), 
            function(object, action, name, name.1, name.2, ...) {
              args <- list(...)
              FactorVertices <- NULL
              FactorEdges <- NULL
              if (!is.null(args$Arguments$ArgBlocks)) 
                warning("Interface for Block-recursive models not implemented!!!")
              f <- function(type)
                if (is.null(type)) 
                  ""
                else paste("(", type, ")")
              if (action == "dropEdge") {
                ## message(paste("Should return an object with the edge from", 
                ## name.1, f(args$from.type), "to", name.2, f(args$to.type), 
                ## "deleted from the argument object"))
                
                new.object <- dropEdge(object,name.1,name.2)
              
              }
              else if (action == "addEdge") {
                ## message(paste("Should return an object with the edge from", 
                ## name.1, f(args$from.type), "to", name.2, f(args$to.type), 
                ## "added to the argument object"))

                new.object <- addEdge(object,name.1,name.2)
                
              }
              else if (action == "dropVertex") {
                ##  message(paste("Should return an object with the vertex", 
                ##  name, f(args$type), "deleted from the argument object"))
                if (!is.null(args$Arguments) &&
                    (args$index > 0)
                    && !is.null(args$Arguments$ArgFactorVertices) && 
                    !is.null(args$Arguments$ArgVertices)) {
                  x <- (args$Arguments$ArgFactorVertices)
                  factors <- lapply(x, function(i) i@vertex.indices)
                  types <- lapply(x, function(i) class(i))
                  factors <- lapply(factors, function(x) x[x != 
                                                           args$index])
                  if (!(is.null(factors))) {
                    result <- returnFactorVerticesAndEdges(args$Arguments$ArgVertices, 
                                                           factors, types)
                    FactorVertices <- result$FactorVertices
                    FactorEdges <- result$FactorEdges
                  }
                }
                
                new.object <- dropVertex(object,name)
              }
              else if (action == "addVertex") {
                ##  message(paste("Should return an object with the vertex", 
                ## name, f(args$type), args$index, "added to the argument object"))
                new.object <- addVertex(object,name)
              
                ## new.object <- subModifyModel(object, action = "add.interactions", 
                ##  modification = name, ...)
              }
              result <- list(object = new.object, FactorVertices = FactorVertices, 
                             FactorEdges = FactorEdges)
              return(result)
            })
  
  if (!isGeneric("testEdge")) {
    if (is.function("testEdge")) 
      fun <- testEdge
    else fun <- function(object, action, name.1, name.2, 
                         ...) standardGeneric("testEdge")
    setGeneric("testEdge", fun)
  }
  setMethod("testEdge", signature(object = "hllm"), 
            function(object, action, name.1, name.2, ...) {
              from.type <- args$from.type
              to.type <- args$to.type
              f <- function(type) if (is.null(type)) 
                ""
              else paste("(", type, ")")
              if (!is.null(args$Arguments$ArgBlocks) || (!is.null(args$Arguments$oriented) && 
                                                         args$Arguments$oriented)) {
                message <- paste("Test of the edge from", name.1, 
                                 "to", name.2, " is not implemented for causal models!!!")
                message(message)
                warning(message)
              }

              if (!extends(class(object),"hllmengine")) 
                object <- fit(object)
              object.small <- dropEdge(object,name.1,name.2)
              
              fit.big <- getFit(object)
              fit.small <- getFit(object.small)
              
              dev.big <- deviance(fit.big)
              df.big  <- fit.big$df
              
              dev.small <- deviance(fit.small)
              df.small  <- fit.small$df
              
              dev.diff <- -(dev.big - dev.small)
              df.diff  <- -(df.big  - df.small)

              return(new("hllmTestClass",
                  deviance = dev.diff,
                  df = df.diff,
                  p = 1-pchisq(dev.diff,df.diff)
                  ))
          })


  ## ####################################################################
  ## end Interface   gRbase <-> dynamicGraph
  ## ####################################################################
}
## Variable types
validVarTypes <- function() c("Discrete", "Ordinal", "Continuous")

## may be extended later
##oldtypes <- validVarTypes()
##validVarTypes <- function() {
##  c(oldtypes,"MyVarType")
##}
.Load.gRbase <- function() {
  .Load.gRbase.general()
  .Load.gRbase.hllm()
  .Load.gRbase.hllmfit()
  .Load.gRbase.hllmodify()
  .Load.dynamicgraph()
}



.First.lib <- function(lib, pkg)
{
  if((R.version$major == 1) && (as.numeric(R.version$minor) < 9))
    packageDescription <- package.description
  
  cat("\n")
  cat("-------------------------------------------------------------\n")
  cat(packageDescription("gRbase", lib = lib, field="Title"))
  cat("\n")
  ver  <- packageDescription("gRbase", lib = lib, field="Version")
  maint<- packageDescription("gRbase", lib = lib, field="Maintainer")
  autho<- packageDescription("gRbase", lib = lib, field="Author")
  descr<- packageDescription("gRbase", lib = lib, field="Description")
  built<- packageDescription("gRbase", lib = lib, field="Built")
  URL  <- packageDescription("gRbase", lib = lib, field="URL")
  cat(descr,"\n")
  cat(paste("gRbase, version", ver,  "is now loaded\n"))
  cat("Authors:",autho,"\n")
  cat("Maintained by",maint,"\n")
  cat("Webpage:",URL,"\n")
  cat("\nBuilt:",built,"\n")
  cat("-------------------------------------------------------------\n")

  require(methods)
  .Load.gRbase()
  
  return(invisible(0))
}

.onAttach <- function (lib, pkg) 
{
    require(methods)
    .Load.gRbase()
  }

.onLoad <- function (lib, pkg) 
{
    require(methods)
    .Load.gRbase()
}


.Last.lib <- function(lib) {
  cat("Thank you for using gRbase\n")
  return(invisible(0))
}
