.packageName <- "graph"

get.circle <- function(points,i,j,r,digits=NULL)
{
	if(is.null(digits)){
		a <- points[i,1]
		b <- points[i,2]
		cc <- points[j,1]
		d <- points[j,2]
	}
	else{
		a <- round(points[i,1],digits=digits)
		b <- round(points[i,2],digits=digits)
		cc <- round(points[j,1],digits=digits)
		d <- round(points[j,2],digits=digits)
	}
	if(a==cc){
	   if(abs(b-d)<=(2*r)){
			 Q <- sqrt(4*r^2-(b-d)^2)
			 x1 <- a-Q/2
			 y1 <- (b+d)/2
			 x2 <- a+Q/2
			 y2 <- (b+d)/2
		  }
		  else return(NULL)
	}
	else if(sqrt((a-cc)^2+(b-d)^2)<=(2*r)){
	  Q <- (a^2+b^2-cc^2-d^2)/(2*(a-cc))
	  R <- (d-b)/(a-cc)
	  S <- Q-a
	  C <- S^2+b^2-r^2
	  B <- 2*S*R-2*b
	  A <- R^2+1
	  Z <- B^2-4*A*C
	  if(abs(Z)<10^(-10)) Z <- 0
	  if(Z>=0){
		 y1 <- (-B+sqrt(Z))/(2*A)
		 x1 <- Q+R*y1
		 y2 <- (-B-sqrt(Z))/(2*A)
		 x2 <- Q+R*y2
	  }
	  else return(NULL)
   }
   else return(NULL)
	rbind(c(x1,y1),c(x2,y2))
}

alpha.hull <- function(points,alpha=NULL,beta=NULL,T=0,digits=NULL)
{
	if(is.null(alpha)){
	   if(is.null(beta)) stop("one of alpha or beta must be given")
		else if(beta<1) stop("beta must be 1 or greater")
	}
   n <- nrow(points)
   d <- ncol(points)
   if(d!=2) stop("alpha hull code only implemented for points in the plane")
   if(is.null(beta)){
		if(is.null(alpha)){
		   m <- mst(points)
			alpha <- sqrt(m$sum/n)
		}
		r <- abs(alpha)
	}
   Adj <- matrix(0,nrow=n,ncol=n)
   if(alpha==0 && is.null(beta)){
      z <- chull(points)
	  n <- length(z)
	  Adj[cbind(z,c(z[2:length(z)],z[1]))] <- 1
   }
   else{
	  for(i in 1:(n-1)){
		 for(j in (i+1):n){
			indices <- setdiff(1:n,c(i,j))
			if(!is.null(beta)){
			   r <- max(as.vector(pdist(points[i,],points[j,],d=2,p=2)*beta)/2,T)
			}
			circs <- get.circle(points,i,j,r,digits)
			if(!is.null(circs)){
			   z <- round(pdist(circs[1,],points[indices,],d=2,p=2),10)
			   if(alpha>0 || !is.null(beta)){
				  if(sum(z<r)==0){
					 Adj[i,j] <- 1
				  }
			   }
			   else{
				  if(sum(z<=r)==length(z)){
					 Adj[i,j] <- 1
				  }
			   }
			   if(Adj[i,j]==0){
				  z <- pdist(circs[2,],points[indices,],d=2,p=2)
				  if(alpha>0 || !is.null(beta)){
					 if(sum(z<r)==0){
						Adj[i,j] <- 1
					 }
				  }
				  else{
					 if(sum(z<=r)==length(z)){
						Adj[i,j] <- 1
					 }
				  }
			   }
			}
		 }
	  }
   }
   g <- as.graph(Adj)
	g$alpha <- alpha
	g$beta <- beta
   g$plot.mode <- "COORDS"
   g$coord <- points
   g
}
beta.skeleton <- function(x,beta=1,p=2,tiny=F)
{
   dx <- pdist(x,p=p,d=ncol(x))
   n <- nrow(x)
   A <- matrix(0,nrow=n,ncol=n)
   if(is.vector(x)) x <- matrix(x,ncol=1)
   for(i in 1:n){
	  for(j in setdiff(1:n,i)){
		 cent1 <- x[i,]*(1-beta/2) + x[j,]*beta/2
		 cent2 <- x[j,]*(1-beta/2) + x[i,]*beta/2
		 r <- beta*dx[i,j]/2
		 d1 <- pdist(x,cent1,d=ncol(x))
		 d2 <- pdist(x,cent2,d=ncol(x))
		 a <- d1<r
		 a[i] <- 0
		 a[j] <- 0
		 b <- d2<r
		 b[i] <- 0
		 b[j] <- 0

		 if(!any(a&b)){
			 A[i,j] <- 1
			 A[j,i] <- 1
		 }
	  }
   }
	if(tiny){
		g <- as.graph(A)
	}
	else {
		g <- as.graph(A)
		g$coord <- x
		g$plot.mode <- "COORDS"
		g$p <- p
		g$beta <- beta
		class(g) <- c("Beta","graph")
	}
   g
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

biggestplaces <- function (x) 
{
	which(x %in% max(x))
}

smallestplaces <- function (x) 
{
	which(x %in% min(x))
}

biggestplace <- function(x)
{
   biggestplaces(x)[1]
}

smallestplace <- function(x)
{
   smallestplaces(x)[1]
}
cccd <- function(x=NULL,y=NULL,dxx=NULL,dyx=NULL,p=2,tiny=F,d)
{
   if(is.null(dxx) | is.null(dyx)){
		if(is.null(x) | is.null(y)) stop("either x,y or dxx,dyx must be given")
		if(missing(d)){
			if(is.matrix(x)){
			   d <- ncol(x)
			}
			else if(is.matrix(y)){
			   d <- ncol(y)
			}
			else d <- 1
		}
      dyx <- pdist(y,x,p=p,d=d)
      dxx <- pdist(x,p=p,d=d)
   }
   R <- apply(dyx,2,min)
	if(tiny) {
		out <- as.digraph(matrix(as.integer(dxx<R),length(R)))
	}
	else{
		out <- as.digraph(matrix(as.integer(dxx<R),length(R)))
		out$R <- R
		out$p <- p
		if(d==1 & !is.null(x) & !is.null(y)){
		   x <- cbind(x,rep(0,length(x)))
		   y <- cbind(y,rep(0,length(x)))
		}
		out$x <- x
		out$y <- y
		class(out) <- c("CCCD","graph")
	}
	out
}

cccd.rw <- function(x=NULL,y=NULL,dxx=NULL,dyx=NULL,p=2,tiny=F,d=2,m=1)
{
   if(is.null(dxx) | is.null(dyx)){
		if(is.null(x) | is.null(y)) stop("either x,y or dxx,dyx must be given")
		if(is.matrix(x)){
			dim <- ncol(x)
		}
		else if(is.matrix(y)){
			dim <- ncol(y)
		}
		else dim <- 1
      dyx <- pdist(y,x,p=p,d=dim)
      dxx <- pdist(x,p=p,d=dim)
   }
   R <- rep(0,nrow(dxx))
	nx <- nrow(dxx)
	ny <- nrow(dyx)
	for(i in 1:nx){
		o <- order(c(dxx[i,],dyx[,i]))
	   rw <- cumsum(c(rep(1/nx,nx),rep(-1/ny,ny))[o])
		r <- sort(dxx[i,])
		R[i] <- r[biggestplace(rw[1:nx]-m*r^d)]
	}
	if(tiny) {
		out <- as.digraph(matrix(as.integer(dxx<R),length(R)))
	}
	else{
		out <- as.digraph(matrix(as.integer(dxx<R),length(R)))
		out$R <- R
		out$p <- p
		if(dim==1 & !is.null(x) & !is.null(y)){
		   x <- cbind(x,rep(0,length(x)))
		   y <- cbind(y,rep(0,length(x)))
		}
		out$x <- x
		out$y <- y
		class(out) <- c("CCCD","graph")
	}
	out
}

plot.CCCD <- function(x,plot.circles=F,dominate.only=F,D=NULL,
                     vertex.pch=19,vertex.col=1,vertex.cex=1,
							xpch=vertex.pch,xcex=vertex.cex,xcol=vertex.col,
							ypch=vertex.pch+1,ycex=vertex.cex,ycol=vertex.col+1,
							circle.col=1, usecoords=T,use.circle.radii=F,
							square=F,
							xlim,ylim,...)
{
	g <- x
	if(usecoords){
	   xcoords <- g$x
	   ycoords <- g$y
		if(missing(xlim)){
			xlim <- range(c(xcoords[,1],ycoords[,1]))
		}
		if(missing(ylim)){
			ylim <- range(c(xcoords[,2],ycoords[,2]))
		}
	}
	else{
	   xcoords <- NULL
	   ycoords <- NULL
	}
	if(dominate.only){
	   if(is.null(D)) D <- dominate(g)
		g$A <- g$A[D,D]
		g$R <- g$R[D]
		if(!is.null(xcoords)){
			x <- xcoords
			xcoords <- xcoords[D,]
		}
	}
	if(use.circle.radii){
	    r <- g$R
		 xlim <- range(c(xlim,xcoords[,1]+r,xcoords[,1]-r))
		 ylim <- range(c(ylim,xcoords[,2]+r,xcoords[,2]-r))
	}
	if(square){
	   xlim <- range(c(xlim,ylim))
	   ylim <- xlim
	}
	if(is.null(xcoords)){
		if(missing(xlim) | missing(ylim)){
			plot.graph(g,vertex.pch=vertex.pch,
						  vertex.col=vertex.col,vertex.cex=vertex.cex,...)
		}
		else{
			plot.graph(g,vertex.pch=vertex.pch,
						  vertex.col=vertex.col,vertex.cex=vertex.cex,
						  xlim=xlim,ylim=ylim,...)
		}
	}
	else{
		plot.graph(g,xlim=xlim,ylim=ylim,vertex.pch=vertex.pch,
					  vertex.col=vertex.col,vertex.cex=vertex.cex,
					  coord=xcoords,...)
		if(dominate.only){
			points(x,pch=xpch,col=xcol,cex=xcex)
		}
		if(!is.null(ycoords)){
			points(ycoords,pch=ypch,col=ycol,cex=ycex)
		}
		if(plot.circles){
			for(i in 1:nrow(xcoords)){
				circle(g$R[i],xcoords[i,1],xcoords[i,2],border=circle.col)
			}
		}
	}
}

cccd.stream <- function(z,cls,a)
{
	sx <- a$sx
	rx <- a$rx
	sy <- a$sy
	ry <- a$ry
	d <- a$d
	clsx <- a$clsx
	if(cls==clsx){
		dx <- as.vector(pdist(z,sx,d=d))
		dy <- as.vector(pdist(z,sy,d=d))
		if(!any(rx>dx)){
		   sx <- rbind(sx,z)
			rx <- c(rx,min(dy))
		}
		ry <- apply(cbind(ry,dy),1,min)
	}
	else{
		dy <- as.vector(pdist(z,sy,d=d))
		dx <- as.vector(pdist(z,sx,d=d))
		if(!any(ry>dy)){
		   sy <- rbind(sy,z)
			ry <- c(ry,min(dx))
		}
		rx <- apply(cbind(rx,dx),1,min)
	}
	list(sx=sx,rx=rx,sy=sy,ry=ry,clsx=clsx,d=d)
}

cccd.classifier <- function(x,y)
{
	if(missing(y)){
	   if(is.list(x)){
		   y <- x$y
			x <- x$x
			if(is.null(x) | is.null(y))
			   stop("must provide either x and y or a list with attributes x and y")
		}
		else
			stop("must provide either x and y or a list with attributes x and y")
	}
	require(graph)
   Gx <- cccd(x,y)
	Gy <- cccd(y,x)
	Dx <- dominate(Gx)
	Dy <- dominate(Gy)
	list(Rx=Gx$R[Dx],Ry=Gy$R[Dy],Cx=matrix(x[Dx,],ncol=ncol(x)),Cy=matrix(y[Dy,],ncol=ncol(y)))
}

cccd.classifier.rw <- function(x,y,m=1,d=2)
{
	if(missing(y)){
	   if(is.list(x)){
		   y <- x$y
			x <- x$x
			if(is.null(x) | is.null(y))
			   stop("must provide either x and y or a list with attributes x and y")
		}
		else
			stop("must provide either x and y or a list with attributes x and y")
	}
	require(graph)
   Gx <- cccd.rw(x,y,m=m,d=d)
	Gy <- cccd.rw(y,x,m=m,d=d)
	Dx <- dominate(Gx)
	Dy <- dominate(Gy)
	list(Rx=Gx$R[Dx],Ry=Gy$R[Dy],Cx=matrix(x[Dx,],ncol=ncol(x)),Cy=matrix(y[Dy,],ncol=ncol(y)))
}

cccd.classify <- function(data,C)
{
	require(graph)
	dx <- apply(t(t(pdist(data,C$Cx))/C$Rx),1,min)
	dy <- apply(t(t(pdist(data,C$Cy))/C$Ry),1,min)
	dx<dy
}
ccd.nonsequential <- function(data,m=1,tiny=F)
{
   r <- rep(0,nrow(data))
   stats <- rep(0,nrow(data))
	walks <- matrix(0,ncol=nrow(data),nrow=nrow(data))
	fs <- matrix(0,ncol=nrow(data),nrow=nrow(data))
	rx <- matrix(0,ncol=nrow(data),nrow=nrow(data))
	for(i in 1:length(r)){
	   y <- data[i,]
	   d <- as.vector(pdistxy(data,y))
	   od <- sort(d)
	   f <- (1:nrow(data))/nrow(data)
	   dif <- f-m*(od/max(od))^2
	   r[i] <- od[biggestplace(dif)]
	   stats[i] <- max(dif)
	   rx[i,] <- od
	   walks[i,] <- f
	   fs[i,] <- (od/max(od))^2
	}
   n <- nrow(data)
   A <- matrix(0,nrow=n,ncol=n)
   for(i in 1:n){
      A[i,] <- pdist(data[i,],data,d=ncol(data))<r[i] 
   }
	if(tiny){
		out <- as.digraph(A)
	}
	else{
		out <- as.digraph(A)
		out$r <- r
		out$stats <- stats
		out$x <- x
		out$walks <- walks
		out$fs <- fs
		out$m <- m
		out$alpha <- NULL
		out$data <- data
		class(out) <- c("CCD","graph")
	}
	out
}

ccd.sequential <- function(data,m=1,alpha=0.05,tiny=F)
{
   n <- nrow(data)
   r <- rep(0,n)
   ks <- sqrt(qchisq(1-alpha,2)/(4*n))
   stats <- rep(0,n)
	walks <- matrix(0,ncol=n,nrow=n)
	fs <- matrix(0,ncol=n,nrow=n)
	rx <- matrix(0,ncol=n,nrow=n)
	for(i in 1:n){
	   y <- data[i,]
	   d <- as.vector(pdistxy(data,y))
	   od <- sort(d)
	   f <- (1:nrow(data))/nrow(data)
	   dif <- f-m*(od/max(od))^2
      # find first one less than -ks
      a <- match(T,c(0,dif)[1:n]< -ks,nomatch=n)
      r[i] <- od[biggestplace(dif[1:a])]
	   stats[i] <- max(dif[1:a])
	   rx[i,] <- od
	   walks[i,] <- f
	   fs[i,] <- (od/max(od))^2
	}
   A <- matrix(0,nrow=n,ncol=n)
   for(i in 1:n){
      A[i,] <- pdist(data[i,],data,d=ncol(data))<r[i] 
   }
	if(tiny)
		out <- as.digraph(A)
	else{
		out <- as.digraph(A)
		out$r <- r
		out$stats <- stats
		out$x <- x
		out$walks <- walks
		out$fs <- fs
		out$m <- m
		out$alpha <- alpha
		out$data <- data
		class(out) <- c("CCD","graph")
   }
	out
}

ccd <- function(data,m=1,alpha=0.05,sequential=T,tiny=F)
{
   if(sequential)
	   z <- ccd.sequential(data,m,alpha,tiny=tiny)
   else
	   z <- ccd.nonsequential(data,m,tiny=tiny)
   z
}

plot.CCD <- function(x,plot.circles=F,dominate.only=F,D=NULL,
                     vertex.pch=19,vertex.col=1,vertex.cex=1,
							pch=vertex.pch,cex=vertex.cex,col=vertex.col,
							circle.col=1, use.coords=T,
							xlim,ylim,...)
{
	g <- x
	if(length(g)==1){
		plot.graph(g,vertex.pch=vertex.pch,
					  vertex.col=vertex.col,vertex.cex=vertex.cex,
					  coord=g$data,...)
	   return
	}
	if(use.coords){
		data <- g$data
		if(missing(xlim)) xlim <- range(g$data[,1])
		if(missing(ylim)) ylim <- range(g$data[,2])
		if(dominate.only){
			if(is.null(D))
				D <- dominate(g)
			g$A <- g$A[D,D]
			g$data <- g$data[D,]
			g$r <- g$r[D]
		}
		plot.graph(g,xlim=xlim,ylim=ylim,vertex.pch=vertex.pch,
					  vertex.col=vertex.col,vertex.cex=vertex.cex,
					  coord=g$data,...)
		if(dominate.only){
			points(data,pch=pch,col=col,cex=cex)
		}
		if(plot.circles){
			for(i in 1:nrow(g$data)){
				circle(g$r[i],g$data[i,1],g$data[i,2],border=circle.col)
			}
		}
	}
	else 
	   plot.graph(g,vertex.pch=vertex.pch,vertex.col=vertex.col,
		           vertex.cex=vertex.cex,...)
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

# draws a circle
circle <- function(r,x,y,len=100,border=1,col)
{
    z <- seq(0,2*pi,length=len)
	z1 <- r*cos(z)+x
	z2 <- r*sin(z)+y
	if(missing(col))
	   polygon(z1,z2,border=border)
	else
	   polygon(z1,z2,border=border,col=col)
}


vote <- function(x)
{
   cls <- unique(x)
	best <- 0
	k <- 0
	for(class in cls){
	   a <- sum(x==class)
	   if(a>best){
		   best <- sum(cls==class)
			k <- class
		}
		else if(a==best){
		   k <- c(k,class)
		}
	}
   k
}

prune <- function(x,classes,prox="Gabriel",ignore.ties=T,...)
{
	MODES <- c("Gabriel","Relative Neighborhood","k-Nearest Neighbor","Minimum Spanning Tree")
	tmp <- charmatch(prox,MODES)
	if(is.null(tmp)){
		 stop("invalid proximity graph or proximity graph not recognized")
	}
	else if(is.na(tmp)){
		 stop("invalid proximity graph or proximity graph not recognized")
	}
	else if(tmp==0){
		 stop("ambiguous proximity graph: retry with more characters")
	}
	else if(tmp<1 || tmp>length(MODES)){
		 stop("invalid proximity graph or proximity graph not recognized")
	}
	mode <- MODES[tmp]
	if(mode=="Gabriel")
		g <- gg(x,...)
	else if(mode=="Relative Neighborhood")
		g <- rng(x,...)
	else if(mode=="k-Nearest Neighbor")
		g <- nng(x,...)
	else if(mode=="Minimum Spanning Tree")
		g <- mst(x,...)
	n <- graph.order(g)
	v <- NULL
	for(i in 1:n){
	   a <- neighborhood(g,v=i,open=T)
		w <- vote(classes[a]) 
		if(ignore.ties){
			if(any(classes[i] %in% w)) v <- c(v,i)
		}
		else {
			if(length(w)==1)
				if(classes[i] == w) v <- c(v,i)
		}
	}
	list(x=x[v,],v=v,graph=g)
}
clique.number <- function(g,U,size=0,m=0)
{
	if(is.digraph(g)) g <- bidirected(g)
   n <- nrow(g$A)
   if(missing(U)){
      od <- degrees(g,cmode="outdegree")
      U <- (1:n)[order(od)]
   }
   if(length(U)==0){
      return(list(max(m,size)))
   }
   while(length(U)>0){
      if((size+length(U))<=m){
         return(m)
      }
      i <- U[1]
      U <- setdiff(U,i)
      N <- (1:n)[neighborhood(g,i,open=T)]
      a <- Recall(g,intersect(U,N),size+1,m)
      m <- max(m,a[[1]])
   }
   m
}

# find all the cliques (brute force)
cliques <- function(g,m){
   if(is.digraph(g)) g <- bidirected(g)
   if(missing(m)){
      m <- clique.number(g)
   }
   od <- degrees(g,cmode="outdegree")+1
   n <- nrow(g$A)
   V <- (1:n)[od>=m]
   if(m==1) return(matrix(V,ncol=1))
   j <- 1
   S <- NULL
   arg <- paste("c(V[i1]",sep="")
   txt <- paste("for(i1 in 1:",length(V)-(m-1),")",sep="")
   for(i in (m-2):0){
      txt <- paste(txt," for(i",j+1," in (i",j,"+1):",length(V)-i,")",sep="")
      arg <- paste(arg,",V[","i",j+1,"]",sep="")
      j <- j+1
   }
   arg <- paste(arg,")",sep="")
   txt <- paste(txt,
      "{",
      "if(is.complete(g$A[",arg,",",arg,"])){S <- rbind(S,",arg,")}",
      "}",
      sep="")
   eval(parse(text=txt))
   S
}

# find one clique (brute force)
clique <- function(g,m){
   if(is.digraph(g)) g <- bidirected(g)
   if(missing(m)){
      m <- clique.number(g)
   }
   else if(m>clique.number(g)) stop("m must be smaller than the clique number")
   od <- degrees(g,cmode="outdegree")+1
   n <- nrow(g$A)
   V <- (1:n)[od>=m]
   if(m==1) return(V[1])
   j <- 1
   S <- NULL
   arg <- paste("c(V[i1]",sep="")
   fnd <- F
   txt <- paste("for(i1 in 1:",length(V)-(m-1),") if(!fnd) ",sep="")
   for(i in (m-2):0){
      txt <- paste(txt," for(i",j+1," in (i",j,"+1):",length(V)-i,")"," if(!fnd) ",sep="")
      arg <- paste(arg,",V[","i",j+1,"]",sep="")
      j <- j+1
   }
   arg <- paste(arg,")",sep="")
   txt <- paste(txt,
      "{",
      "if(is.complete(g$A[",arg,",",arg,"])){S <- rbind(S,",arg,"); fnd<-T}",
      "}",
      sep="")
   eval(parse(text=txt))
   S
}

cliques.dist <- function(g)
{
   m <- clique.number(g)
   dist <- rep(0,m)
   for(i in 1:m){
      dist[i] <- nrow(cliques(g,i))
   }
   dist
}


declique <- function(G)
{
	A <- adjacency.matrix(G)
	n <- nrow(A)
   if(is.digraph(G)){
		G <- bidirected(G)
	}
	q <- clique.number(G)
	if(q > 2){
		if(q==n) return(as.graph(matrix(0,ncol=1,nrow=1)))
	   a <- clique(G)
		inds <- setdiff(1:n,a)
		B <- matrix(0,nrow=length(inds)+1,ncol=length(inds)+1)
		B[1:length(inds),1:length(inds)] <- A[inds,inds]
		for(i in 1:length(inds)){
		   B[length(inds)+1,i] <- any(A[a,inds[i]])
		   B[i,length(inds)+1] <- any(A[inds[i],a])
		}
	}
	else{
	  B <- A
	}
	G$A <- B
	G$clique <- a
	G
}


clustering.coef <- function(g,k=1,type=1)
{
	if(type==1){
		out <- mean(local.density(g,k))
	}
	else if(type==1){
		out <- 0
		s <- 0
		for(i in 1:graph.order(g)){
			N <- hop.disk(g,v=i,hops=k)
			if(length(N)>1){
				a <- choose(length(N),2)
				out <- out + graph.size(subgraph(g,N))
				s <- s+a
			}
		}
		out <- out/s
	}
	else{
	   stop("type must be 1 or 2")
	}
	out
}
collapse <- function(G,vertices=clique(G))
{
	A <- adjacency.matrix(G)
	n <- nrow(A)
	q <- length(vertices)
	if(q > 1){
		if(q==n) return(as.graph(matrix(0,ncol=1,nrow=1)))
		inds <- setdiff(1:n,vertices)
		B <- matrix(0,nrow=length(inds)+1,ncol=length(inds)+1)
		B[1:length(inds),1:length(inds)] <- A[inds,inds]
		for(i in 1:length(inds)){
		   B[length(inds)+1,i] <- any(A[vertices,inds[i]])
		   B[i,length(inds)+1] <- any(A[inds[i],vertices])
		}
	}
	else{
	   B <- A
	}
	G$A <- B
	G$collapsed <- vertices
	G
}

greedy.color <- function(g)
{
   n <- graph.order(g)
   vc <- rep(0,n)
   vc[1] <- 1
   for(i in 2:n){
      x <- neighborhood(g,i,as.indices=T)
	  if(is.null(x)) vc[i] <- 1
	  else vc[i] <- min((1:n)[-vc[x]])
   }
   vc
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

complexity <- function(x=NULL,dx=NULL,classes,p=2,tiny=F)
{
	if(missing(classes)) stop("classes must be given.")
	if(is.null(dx)) {
	   if(is.null(x))
			stop("x or dx must be given.")
	   dx <- pdist(x,p=p)
	}
   G <- mst(distmat=dx)
	n <- nrow(dx)
	edges <- G$edges
	com <- 0
	A <- matrix(0,nrow=nrow(dx),ncol=ncol(dx))
	for(i in 1:G$nedges){
	   if(classes[edges[i,1]] != classes[edges[i,2]]){
			com <- com+1
			A[edges[i,1],edges[i,2]] <- 1
		}
	}
	if(tiny){
		out <- as.graph(A)
		out$complexity <- com/n
	}
	else{
		out <- as.graph(A)
		out$complexity <- com/n
		out$coord <- x
		out$plot.mode <- "COORDS"
		out$dx <- dx
		out$classes <- classes
		out$p <- p
	}
	class(out) <- c("CG","graph")
	out
}

is.connected <- function(g,mode="weak")
{
	if(mode=="weak"){
	   g <- symmetrize.graph(g)
	}
   A <- adjacency.matrix(g)
	n <- nrow(A)
	marked <- rep(F,n)
	q <- 1
	marked[q] <- T
	all <- 1:n
	while(length(q)>0){
	   w <- q[1]
		q <- setdiff(q,w)
		N <- all[neighborhood(g,w)]
		v <- N[!marked[N]]
		marked[v] <- T
		q <- c(q,v)
	}
	!any(!marked)
}

graph.components <- function(g)
{
	g <- symmetrize.graph(g)
	covered <- adjacency.matrix(g,augmented=T)
   n <- nrow(covered)
   comp <- rep(0,n)
   nc <- apply(covered,1,sum)
   curcomp <- 1
   while(any(comp==0)){
	  i <- biggestplace(nc)
	  osum <- n
      comp[covered[i,]>0] <- curcomp
	  while(osum != sum(comp==curcomp)){
		 cur <- comp==curcomp
		 osum <- sum(cur)
		 if(osum==1){
			v <- covered[cur,]
		 }
		 else{
			v <- apply(covered[cur,],2,max) 
		 }
		 comp[v>0] <- curcomp
	  }
	  if(osum!=n) {
		 nc[comp>0] <- 0
		 curcomp <- curcomp+1
	  }
   }
   comp
}

number.components <- function(g)
{
   max(graph.components(g))
}

component.sequence <- function(g)
{
   sort(table(graph.components(g)))
}

graph.component <- function(g,component=0)
{
   comp <- graph.components(g)
	if(component==0){
	   h <- table(comp)
	   component <- which(h == max(h))
	}
	if(is.digraph(g)){
	   g <- as.digraph(adjacency.matrix(g)[comp==component,comp==component])
	}
	else{
	   g <- as.graph(adjacency.matrix(g)[comp==component,comp==component])
	}
	g
}
delaunay <- function(z)
{
   require(deldir)
   del <- deldir(z[,1],z[,2])
   d <- del$delsgs[,5:6]

   A <- matrix(0,nrow=nrow(z),ncol=nrow(z))
   for(i in 1:nrow(d))
	  A[d[i,1],d[i,2]] <- 1

   g <- as.graph(A)

   g$plot.mode <- "COORDS"

   g$coord <- z

   g
}

nn.edit <- function(x,classes)
{
   require(deldir)
   g <- delaunay(x)
   marks <- rep(F,graph.order(g))
   for(i in 1:length(marks)){
      q <- neighborhood(g,i,open=F,as.indices=T)
	  cls <- classes[q]
	  if(length(unique(cls))>1) marks[i] <- F
   }
   list(x=x[marks,],classes=classes[marks],marks=marks)
}
dominate.greedy <- function(g,weight=NULL)
{
   od <- degrees(g,cmode="outdegree")+1
   S <- NULL
   A <- g$A
   diag(A) <- 1
   n <- nrow(g$A)
   covered <- rep(0,n)
   while(sum(covered)<n){
      i <- match(max(od),od)
      if(!is.null(weight)){
         qq <- weight
         qq[od!=od[i]] <- -Inf
         i <- match(max(weight[od==od[i]]),qq)
      }
      covered[A[i,]>0] <- 1
      S <- c(S,i)
      A[,covered>0] <- 0
		h <- as.digraph(A)
      od <- degrees(h,cmode="outdegree")+1-covered
   }
   S
}

dominate.byR <- function(g,R,eliminate.covered=F)
{
   n <- nrow(g$A)
   S <- NULL
   covered <- rep(F,n)
	A <- g$A
   diag(A) <- 1
   while(any(!covered)){
      v <- match(max(R),R)
      S <- c(S,v)
      cov <- A[v,]==1
      covered <- covered | cov
      R[v] <- -Inf
      if(eliminate.covered)
         R[cov] <- -Inf
   }
   S
}

dominate.random <- function(g,numits=100,prob=0.5)
{
   S <- dominate.greedy(g)
   gamma <- length(S)
   if(gamma <= 2) return(S)
	A <- g$A
   n <- nrow(A)
   zeros <- indegree0(g)
   W <- (1:n)[zeros==F]
   Z <- (1:n)[zeros]
   if(sum(neighborhood(A,Z))==n) return(Z)
   if(sum(zeros)<n){
      for(i in 1:numits){
         q <- gamma-length(Z)-1
         if(q==1)
            X <- c(Z,sample(W,1))
         else
            X <- c(Z,sample(W,rbinom(1,q-1,prob)+1))
         if(sum(neighborhood(A,X))==n){
            S <- X
            gamma <- length(X)
         }
      }
   }
   S
}

dominate.systematic <- function(g,S,verbose=T)
{
	x <- degrees(g,cmode="outdegree")
	if(max(x) == nrow(g$A)) return(match(max(x),x))
   if(missing(S)) S <- dominate.greedy(g)
   gamma <- length(S)
   if(verbose){
      cat("Starting with gamma =",gamma,"\n")
   }
   if(gamma <= 2) return(S)
	A <- g$A
   n <- nrow(A)
   zeros <- indegree0(g)
   Z <- (1:n)[zeros]
   W <- (1:n)[zeros==F]
   if(gamma-1-length(Z)==0){
      if(sum(neighborhood(g,Z))==n){
         return(Z)
      }
      else{
         return(S)
      }
   }
   if(sum(zeros)<n){
      numchoices <- choose(length(W),gamma-1-length(Z))
      if(verbose){
         cat("This may require","(",length(W),",",gamma-1-length(Z),") =",
             numchoices,"steps\n")
      }
      m <- gamma-1-length(Z)
      txt <- paste("for(i1 in 1:",length(W)-(m-1),")",sep="")
      arg <- paste("W[i1]",sep="")
      j <- 1
      fnd <- 0
      if(m>1){
         for(i in (m-2):0){
            txt <- paste(txt," if(!fnd) for(i",j+1," in (i",j,"+1):",length(W)-i,")",sep="")
            arg <- paste(arg,",W[","i",j+1,"]",sep="")
            j <- j+1
         }
      }
      txt <- paste(txt,
         "{",
         "if(!fnd){",
         " S1 <- c(Z,",arg,");",
         "if(sum(neighborhood(g,S1))==",n,"){fnd <- 1; break }",
         "}",
         "}",
      sep="")
      eval(parse(text=txt))
      if(fnd){
         S <- S1
         if(verbose){
            cat("Recursing:",S,"\n")
         }
         Q <- Recall(g,S)
			if(length(Q)<length(S)) S <- Q
      }
   }
   S
}

dominate.all <- function(g,gamma,verbose=T)
{
	A <- g$A
   n <- nrow(A)
   if(missing(gamma)) {
      gamma <- length(dominate.greedy(g))
      if(verbose){
         cat("Setting gamma to the greedy value of",gamma,"\n")
      }
   }
   if(gamma==1){
      a <- degrees(g,cmode="outdegree")
      if(any(a==(n-1))){
         return((1:n)[a==(n-1)])
      }
      else{
         return(NULL)
      }
   }
   zeros <- indegree0(g)
   Z <- (1:n)[zeros]
   W <- (1:n)[zeros==F]
   if(length(Z)==gamma) return(Z)
   if(length(Z)>gamma) return(NULL)
   if(sum(zeros)<n){
      numchoices <- choose(length(W),gamma-length(Z))
      if(verbose){
         cat("This will require","(",length(W),",",gamma-length(Z),") =",
             numchoices,"steps\n")
      }
      m <- gamma-length(Z)
      txt <- paste("for(i1 in 1:",length(W)-(m-1),")",sep="")
      arg <- paste("W[i1]",sep="")
      j <- 1
      doms <- NULL
      if(m>1){
         for(i in (m-2):0){
            txt <- paste(txt," for(i",j+1," in (i",j,"+1):",length(W)-i,")",sep="")
            arg <- paste(arg,",W[","i",j+1,"]",sep="")
            j <- j+1
         }
      }
      txt <- paste(txt,
         "{",
         " S1 <- c(Z,",arg,");",
         "if(sum(neighborhood(g,S1))==",n,"){doms <- rbind(doms,S1)}",
         "}",
         sep="")
         eval(parse(text=txt))
   }
   dimnames(doms) <- NULL
   doms
}

dominate.dist <- function(g,mingamma=1,maxgamma=nrow(g))
{
   if(mingamma>maxgamma) stop("min and max gammas must be chosen appropriately")
   n <- maxgamma
   doms <- rep(0,maxgamma-mingamma+1)
   for(i in mingamma:n){
      a <- dominate.all(g,i,verbose=F)
      if(!is.null(a)) {
         if(is.vector(a)) doms[i-mingamma+1] <- length(a)
         else{
            doms[i-mingamma+1] <- nrow(a)
         }
      }
   }
   doms
}

dominate <- function(g,method="greedy",tiebreak=NULL,
                     verbose=F,...)
{
   METHODS = c("greedy","systematic","random","byR","all")
   method <- pmatch(method,METHODS)
   if(is.na(method)){
      stop("invalid method")
   }
   if(method==1) S <- dominate.greedy(g,tiebreak)
   else if(method==2) S <- dominate.systematic(g,verbose=verbose)
   else if(method==3) S <- dominate.random(g,...)
   else if(method==4) S <- dominate.byR(g,tiebreak)
   else if(method==5){
      S <- dominate.random(g,...)
      S <- dominate.systematic(g,S,verbose=verbose)
   }
   S
}

dominate.random.sample <- function(g,method="greedy")
{
   n <- graph.order(g)
   METHODS = c("greedy","degree")
   method <- pmatch(method,METHODS)
   if(is.na(method)){
      stop("invalid method")
   }
   if(method==1) {
		S <- dominate.greedy(g,runif(n))
	}
	if(method==2){
		S <- NULL
		A <- adjacency.matrix(g)
		diag(A) <- 1
		covered <- rep(0,n)
		od <- degrees(g,cmode="outdegree")+1
		while(sum(covered)<n){
			i <- sample(n,1,prob=od/sum(od))
			covered[A[i,]>0] <- 1
			S <- c(S,i)
			A[,covered>0] <- 0
			h <- as.digraph(A)
			od <- degrees(h,cmode="outdegree")+1-covered
		}
	}
	S
}

dominate.all.greedy <- function(g)
{
   .dominate.all.greedy(g)
	out <- .tmp.greedy.xxx
	a <- list(0)
	k <- 1
	for(i in 1:nrow(out)) {
			new <- T
			S <- sort(out[i,])
			for(j in 1:length(a)){
				if(length(a[[j]]) == length(S)){
					if(all(S==sort(a[[j]]))){
						new <- F
						break
					}
				}
			}
			if(new){
				a[[k]] <- which(out[i,])
				k <- k+1
			}
	}
	a
}

.dominate.all.greedy <- function(g,
	n=graph.order(g),
	covered=rep(0,n),
	S=rep(F,n),
	A=adjacency.matrix(g),
	od=NULL)
{
	if(sum(covered) == n) {
		.tmp.greedy.xxx <<- rbind(.tmp.greedy.xxx,S)
		return(S)
	}
	if(is.null(od)){
		od <- degrees(g,cmode="outdegree")+1
		diag(A) <- 1
		.tmp.greedy.xxx <<- NULL
	}
	if(any(od==1)){
		S[which(od==1)] <- T
		covered[S] <- 1
	}
	a <- which(max(od)==od)
	for(i in 1:length(a)){
		nS <- S
		nS[a[i]] <- T
		#print(nS)
		ncovered <- covered
		ncovered[A[a[i],]>0] <- 1
		nA <- A
		nA[,ncovered>0] <- 0
		h <- as.digraph(nA)
		nod <- degrees(h,cmode="outdegree")+1-ncovered
	   #print(Recall(g,n,ncovered,nS,nA,nod))
	   Recall(g,n,ncovered,nS,nA,nod)
	}
}

odominate.all.greedy <- function(g)
{
   od <- degrees(g,cmode="outdegree")+1
	S <- NULL
   A <- adjacency.matrix(g)
   diag(A) <- 1
   n <- graph.order(g)
   covered <- rep(0,n)
	if(any(od==1)){
		S <- which(od==1)
		covered[S] <- 1
	}
	out <- list(list(S=S,A=A,od=od,covered=covered))
	k <- 1
	kk <- 1
	done <- F
   while(!done){
		covered <- out[[k]]$covered
		A <- out[[k]]$A
		S <- out[[k]]$S
		od <- out[[k]]$od
		while(sum(covered)<n){
			a <- which(max(od)==od)
			if(length(a)>1){
				for(j in 2:length(a)){
					new <- T
					SS <- sort(c(S,a[j]))
					for(jj in 1:length(out)){
						if(length(out[[jj]]$S) == length(SS)){
							if(all(SS==sort(out[[jj]]$S))){
								new <- F
								break
							}
						}
					}
					if(new){
						kk <- kk+1
						out[[kk]] <- list(S=S,A=A,covered=covered)
						out[[kk]]$covered[out[[kk]]$A[a[j],]>0] <- 1
						out[[kk]]$S <- c(out[[kk]]$S,a[j])
						out[[kk]]$A[,out[[kk]]$covered>0] <- 0
						h <- as.digraph(out[[kk]]$A)
						out[[kk]]$od <- degrees(h,cmode="outdegree")+1-out[[kk]]$covered
					}
				}
			}
			i <- a[1]
			covered[A[i,]>0] <- 1
			S <- c(S,i)
			A[,covered>0] <- 0
			h <- as.digraph(A)
			od <- degrees(h,cmode="outdegree")+1-covered
		}
		cat(k,kk,"\n")
		out[[k]] <- list(S=S,A=A,covered=covered)
		if(k==length(out)) done<-T
		else k <- k+1
	}
	a <- list(0)
	k <- 1
	for(i in 1:length(out)) {
			new <- T
			S <- sort(out[[i]]$S)
			for(j in 1:length(a)){
				if(length(a[[j]]) == length(S)){
					if(all(S==sort(a[[j]]))){
						new <- F
						break
					}
				}
			}
			if(new){
				a[[k]] <- out[[i]]$S
				k <- k+1
			}
	}
   a
}

mkdot <- function(G,file="/tmp/dg",labels=NULL,width=.1)
{
   write("digraph G {",file=file)
	n <- nrow(G$A)
	write(paste("node[width =",width,"]"),file=file,append=T)
	for(i in 1:n){
		if(is.null(labels))
			write(paste(i,"[label=\"\",shape=circle];"),file=file,append=T)
		else
			write(paste(i,"[label =",labels[i],"];"),file=file,append=T)
	}
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			if(G$A[i,j]){
				if(G$A[j,i]){
					write("edge[dir=both];",file=file,append=T)
					write(paste(i,"->",j,";"),file=file,append=T)
					write("edge[dir=forward];",file=file,append=T)
				}
				else
					write(paste(i,"->",j,";"),file=file,append=T)
			}
			else if(G$A[j,i]){
					write(paste(j,"->",i,";"),file=file,append=T)
			}
		}
	}
	write("}",file=file,append=T)
}

dotlayout <- function(G)
{
	file <- tempfile()
	n <- nrow(G$A)
   mkdot(G,file)
	df <- paste(file,"dot",sep=".")
	system(paste("dot -o",df,file))
	unlink(file)
	a <- scan(df,what="")
	unlink(df)
	coords <- matrix(0,nrow=n,ncol=2)
	i <- 1
	j <- 1
	while(i<length(a) && j <= n){
	   if(a[i]==paste(j)){
			i <- i+1
			while(i < length(a) && length(grep("pos",a[i]))==0) i <- i+1
			b <- strsplit(a[i],"\"")[[1]][2]
			coords[j,] <- as.integer(strsplit(b,",")[[1]])
		   j <- j+1
		}
		i <- i+1
	}
	coords
}
read.erdos <- function(file=paste(R.home(),"/library/graph/data/erdos1graph",sep=""))
{
   a <- readLines(file)
	n <- as.integer(a[[22]])
	A <- matrix(0,nrow=n,ncol=n)
	names <- rep(0,n)
	for(i in 1:n){
	   b <- a[[22+2*(i-1)+1]]
		if(a[[22+2*i]]!=""){
			x <- as.integer(strsplit(a[[22+2*i]]," ")[[1]])
			x <- x[!is.na(x)]
			if(length(x)>0){
				A[i,x] <- 1
				A[x,i] <- 1
			}
		}
		y <- substring(b,16)
		names[i] <- sub(' +$', '', y)
	}
	E <- as.graph(A)
	E$names <- names
	E
}

erdos <- function() read.erdos()

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

gamma1 <- function(x,y,k=5,p=2,mutual=F,tiny=F)
{
	if(missing(x)) stop("x must be given.")
	if(is.matrix(x)){
	   y <- x[,2]
		x <- x[,1]
	}
	else if(is.list(x)){
	   x <- x$x
		y <- x$y
		if(is.null(x) || is.null(y)){
		   stop("Either x must be a matrix, a list with attributes x and y, or a vector, with y given.")
		}
	}
	else{
	   if(missing(y)) stop("Either x must be a matrix, or y must be given.")
	}
   xG <- symmetrize.graph(nng(x,k=k,p=p,mutual=mutual))
   yG <- symmetrize.graph(nng(y,k=k,p=p,mutual=mutual))
	A <- xG$A & yG$A
	if(tiny){
		out <- as.graph(A)
		out$gamma1 <- sum(A)/2
	}
	else{
		out <- as.graph(A)
		out$gamma1 <- sum(A)/2
		out$k <- k
		out$p <- p
		out$mutual <- mutual
		out$coord <- x
		out$plot.mode <- "COORDS"
		out$y <- y
	}
	class(out) <- c("gamma1","graph")
	out
}


# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

gamma2 <- function(x,y,k=5,p=2,mutual=F)
{
	if(missing(x)) stop("x must be given.")
	if(is.matrix(x)){
	   y <- x[,2]
		x <- x[,1]
	}
	else if(is.list(x)){
	   x <- x$x
		y <- x$y
		if(is.null(x) || is.null(y)){
		   stop("Either x must be a matrix, a list with attributes x and y, or a vector, with y given.")
		}
	}
	else{
	   if(missing(y)) stop("Either x must be a matrix, or y must be given.")
	}
   xG <- symmetrize.graph(nng(x,k=k,p=p,mutual=mutual))
	A <- xG$A
	n <- nrow(A)
	g <- 0
	if(is.vector(x)) x <- matrix(x,ncol=1)
	if(is.vector(y)) y <- matrix(y,ncol=1)
	d <- ncol(y)
	for(i in 1:(n-1)){
	   for(j in (i+1):n){
		   if(A[i,j]>0){
			   dd <- pdist(y[i,],y,p=p,d=d)
				g <- g+order(dd)[j]
			}
		}
	}
	g
}


# uses Floyd-Warshal, Jungnickel, p. 87
geodesic.dist <- function(g)
{
   d <- adjacency.matrix(g)
	n <- nrow(d)
	d[d==0] <- n+2
	diag(d) <- 0
	for(k in 1:n){
		for(i in 1:n){
			for(j in 1:n){
				d[i,j] <- min(d[i,j],d[i,k]+d[k,j])
			}
		}
	}
	d
}
gg <- function(x,r=1,p=2,tiny=F,usedeldir=T)
{
	dx <- pdist(x,p=p)
   n <- nrow(dx)
   A <- matrix(0,nrow=n,ncol=n)
   if(is.vector(x)) x <- matrix(x,ncol=1)
   if(usedeldir && p==2 && ncol(x)==2 && (length(.find.package("deldir",quiet=T))>0)){
	  require(deldir)
	  del <- deldir(x[,1],x[,2])
	  for(edge in 1:nrow(del$delsgs)){
	    i <- del$delsgs[edge,5]
	    j <- del$delsgs[edge,6]
		 d1 <- r*dx[i,j]/2
		 d <- pdist((x[i,]+x[j,])/2,x,d=ncol(x),p=p)
		 d[i] <- Inf
		 d[j] <- Inf
		 if(!any(d<d1)){
		    A[i,j] <- 1
		    A[j,i] <- 1
		 }
	  }
   }
   else{
	  for(i in 1:n){
		 for(j in setdiff(1:n,i)){
			d1 <- r*dx[i,j]/2
			d1 <- round(d1,10)
			d <- pdist((x[i,]+x[j,])/2,x,d=ncol(x),p=p)
			d <- round(d,10)
			d[i] <- Inf
			d[j] <- Inf
			if(!any(d<d1)){
				A[i,j] <- 1
				A[j,i] <- 1
			}
		 }
	  }
   }
	if(tiny){
		g <- as.graph(A)
	}
	else {
		g <- as.graph(A)
		g$r <- r
		g$coord <- x
		g$plot.mode <- "COORDS"
		g$p <- p
		class(g) <- c("GG","graph")
	}
   g
}

# petersen graph
petersen <- function(){
	A <- rbind(
				 c(0,1,0,0,1,1,0,0,0,0),
				 c(1,0,1,0,0,0,1,0,0,0),
				 c(0,1,0,1,0,0,0,1,0,0),
				 c(0,0,1,0,1,0,0,0,1,0),
				 c(1,0,0,1,0,0,0,0,0,1),
				 c(1,0,0,0,0,0,0,1,1,0),
				 c(0,1,0,0,0,0,0,0,1,1),
				 c(0,0,1,0,0,1,0,0,0,1),
				 c(0,0,0,1,0,1,1,0,0,0),
				 c(0,0,0,0,1,0,1,1,0,0))
	out <- as.graph(A)
	out$plot.mode <- "mcircles"
	out$m <- 2
	out
}

star <- function(n)
{
   A <- matrix(0,nrow=n,ncol=n)
	for(i in 2:n){
	   A[1,i] <- 1
	   A[i,1] <- 1
	}
	g <- as.graph(A)
	g$plot.mode <- "COORDS"
	g$coord <- cbind(c(0,sin(2 * pi * ((0:(n - 2))/(n-1)))),
						  c(0,cos(2 * pi * ((0:(n - 2))/(n-1)))))
   g
}

wheel <- function(n)
{
	g1 <- star(n)
	g2 <- Cn(n-1,gmode="graph")
   A <- adjacency.matrix(g1)
	A[2:n,2:n] <- adjacency.matrix(g2)
	g <- as.graph(A)
	g$plot.mode <- "COORDS"
	g$coord <- g1$coord
   g
}

fan <- function(n)
{
   g <- wheel(n)
	g$A[2,n] <- 0
	g$A[n,2] <- 0
	g$coord <- cbind(c(0,cos(seq(0+.1,pi-.1,length=n-1))),
						  c(0,sin(seq(0+.1,pi-.1,length=n-1))))
   g
}

# function to return a complet graph
Kn <- function(n){
   g <- matrix(1,nrow=n,ncol=n)
   diag(g) <- 0
	as.graph(g)
}

# a function to return a cycle
Cn <- function(n,gmode="digraph"){
   A <- matrix(0,nrow=n,ncol=n)
   for(i in 1:(n-1)) A[i,i+1] <- 1
   A[n,1] <- 1
	g <- as.graph(A)
   if(gmode=="graph") g <- symmetrize.graph(g)
	g$plot.mode <- "circle"
   g
}
# a function to return a chain
chain <- function(n,gmode="digraph"){
   A <- matrix(0,nrow=n,ncol=n)
	for(i in 1:(n-1)){
	   A[i,i+1] <- 1
		if(gmode=="graph")
			A[i+1,i] <- 1
	}
	g <- as.graph(A)
	g$plot.mode="line"
	g
}

# complete bipartite graph
Knm <- function(n,m){
	A <- matrix(0,nrow=(n+m),ncol=(n+m))
	A[1:n,(n+1):(n+m)] <- 1
	A[(n+1):(n+m),1:n] <- 1
	out <- as.graph(A)
	out$plot.mode <- "bipartite"
	out$parts <- c(rep(1,n),rep(2,m))
	class(out) <- c("graph","bipartite")
	out
}

trivial.graph <- function()
{
	out <- as.graph(matrix(0,1,1))
	out
}

empty.graph <- function(n)
{
   as.graph(matrix(0,nrow=n,ncol=n))
}

tetrahedron <- function()
{
   g <- Kn(4)
   coord <- rbind(
         c(0,1),
         c(-1,0),
         c(1,0),
         c(0,.35))
   g$plot.mode <- "COORDS"
   g$coord <- coord
   g
}

hypercube <- function(d=3)
{
   n <- 2^d
	A <- matrix(0,nrow=n,ncol=n)
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			out <- 0
			x <- i-1
			y <- j-1
			for(k in 1:d){
				z <- floor(x/2)
			   a <- abs(x-2*z)
				x <- z
				z <- floor(y/2)
			   b <- abs(y-2*z)
				y <- z
				out <- out + (a==b)
			}
		   if(out==(d-1)){
			   A[i,j] <-1
			   A[j,i] <-1
			}
		}
	}
	as.graph(A)
}

square <- function()
{
   g <- hypercube(2)
	coords <- rbind(c(0,0),c(0,1),c(1,0),c(1,1))
	g$plot.mode <- "COORDS"
	g$coord <- coords
	g
}

cube <- function()
{
   g <- hypercube(3)
	coords <- rbind(c(0,0),c(0,1),c(1,0),c(1,1),
	                c(-.5,-.5),c(-.5,1.5),c(1.5,-.5),c(1.5,1.5))
	g$plot.mode <- "COORDS"
	g$coord <- coords
	g
}

octahedron <- function()
{
   n <- 6
	A <- matrix(0,nrow=n,ncol=n)
	A[1,c(2,3,4,5)] <- 1
	A[2,c(1,3,4,6)] <- 1
	A[3,c(1,2,5,6)] <- 1
	A[4,c(1,2,5,6)] <- 1
	A[5,c(1,3,4,6)] <- 1
	A[6,c(2,3,4,5)] <- 1
	A[c(2,3,4,5),1] <- 1
	A[c(1,3,4,6),2] <- 1
	A[c(1,2,5,6),3] <- 1
	A[c(1,2,5,6),4] <- 1
	A[c(1,3,4,6),5] <- 1
	A[c(2,3,4,5),6] <- 1
	g <- as.graph(A)
	g$plot.mode <- "COORDS"
	coords <- rbind(c(-1,1),c(1,1),c(0,0),
	                c(0,.75),c(-.25,.5),c(.25,.5))
   g$coord <- coords
	g
}

dodecahedron <- function()
{
	n <- 20
	A <- matrix(0,nrow=n,ncol=n)
	A[1,c(2,5,6)] <- 1
	A[c(2,5,6),1] <- 1
	A[2,c(1,3,8)] <- 1
	A[c(1,3,8),2] <- 1
	A[3,c(2,4,10)] <- 1
	A[c(2,4,10),3] <- 1
	A[4,c(3,5,12)] <- 1
	A[c(3,5,12),4] <- 1
	A[5,c(1,4,14)] <- 1
	A[c(1,4,14),5] <- 1
	A[6,c(1,7,15)] <- 1
	A[c(1,7,15),6] <- 1
	A[7,c(6,8,17)] <- 1
	A[c(6,8,17),7] <- 1
	A[8,c(2,7,9)] <- 1
	A[c(2,7,9),8] <- 1
	A[9,c(8,10,18)] <- 1
	A[c(8,10,18),9] <- 1
	A[10,c(3,9,11)] <- 1
	A[c(3,9,11),10] <- 1
	A[11,c(10,12,19)] <- 1
	A[c(10,12,19),11] <- 1
	A[12,c(4,11,13)] <- 1
	A[c(4,11,13),12] <- 1
	A[13,c(12,14,20)] <- 1
	A[c(12,14,20),13] <- 1
	A[14,c(5,13,15)] <- 1
	A[c(5,13,15),14] <- 1
	A[15,c(6,14,16)] <- 1
	A[c(6,14,16),15] <- 1
	A[16,c(15,17,20)] <- 1
	A[c(15,17,20),16] <- 1
	A[17,c(7,16,18)] <- 1
	A[c(7,16,18),17] <- 1
	A[18,c(9,17,19)] <- 1
	A[c(9,17,19),18] <- 1
	A[19,c(11,18,20)] <- 1
	A[c(11,18,20),19] <- 1
	A[20,c(13,16,19)] <- 1
	A[c(13,16,19),20] <- 1
   coord1 <- cbind(3.5*sin(2 * pi * ((0:4)/5)),
						 3.5*cos(2 * pi * ((0:4)/5)))
   coord2 <- cbind(2*sin(2 * pi * ((0:9)/10)),
						 2*cos(2 * pi * ((0:9)/10)))
   coord3 <- cbind(sin(2 * pi * (((-1):3)/5+.1)),
						 cos(2 * pi * (((-1):3)/5+.1)))
	g <- as.graph(A)
   g$plot.mode="COORDS"
	g$coord=rbind(coord1,coord2,coord3)
	g
}

icosahedron <- function()
{
   n <- 12
	A <- matrix(0,nrow=n,ncol=n)
	A[1,c(2,3,4,8,9)] <- 1
	A[2,c(1,3,4,5,6)] <- 1
	A[3,c(1,2,6,7,8)] <- 1
	A[4,c(1,2,5,9,10)] <- 1
	A[5,c(2,4,6,10,11)] <- 1
	A[6,c(2,3,5,7,11)] <- 1
	A[7,c(3,6,8,11,12)] <- 1
	A[8,c(1,3,7,9,12)] <- 1
	A[9,c(1,4,8,10,12)] <- 1
	A[10,c(4,5,9,11,12)] <- 1
	A[11,c(5,6,7,10,12)] <- 1
	A[12,c(7,8,9,10,11)] <- 1
	A[c(2,3,4,8,9),1] <- 1
	A[c(1,3,4,5,6),2] <- 1
	A[c(1,2,6,7,8),3] <- 1
	A[c(1,2,5,9,10),4] <- 1
	A[c(2,4,6,10,11),5] <- 1
	A[c(2,3,5,7,11),6] <- 1
	A[c(3,6,8,11,12),7] <- 1
	A[c(1,3,7,9,12),8] <- 1
	A[c(1,4,8,10,12),9] <- 1
	A[c(4,5,9,11,12),10] <- 1
	A[c(5,6,7,10,12),11] <- 1
	A[c(7,8,9,10,11),12] <- 1
	g <- as.graph(A)
   g$plot.mode="COORDS"
	coord1 <- rbind(c(-6,6),c(6,6),c(0,-4))
	coord2 <- cbind(2*sin(2 * pi * ((0:(7 - 2))/(7-1))),
						 2*cos(2 * pi * ((0:(7 - 2))/(7-1)))+2)
	coord3 <- rbind(c(0,3),c(1,1.5),c(-1,1.5))
	g$coord=rbind(coord1,coord2,coord3)
	g
}

birkhoff.diamond <- function()
{
   n <- 10
	A <- matrix(0,nrow=n,ncol=n)
	A[1,c(2,6,7)] <- 1
	A[2,c(1,3,7,8)] <- 1
	A[3,c(2,4,8,9)] <- 1
	A[4,c(3,5,9)] <- 1
	A[5,c(4,6,9,10)] <- 1
	A[6,c(1,5,7,10)] <- 1
	A[7,c(1,2,6,8,10)] <- 1
	A[8,c(2,3,7,9)] <- 1
	A[9,c(3,4,5,8,10)] <- 1
	A[10,c(5,6,7,9)] <- 1
	A[c(2,6,7),1] <- 1
	A[c(1,3,7,8),2] <- 1
	A[c(2,4,8,9),3] <- 1
	A[c(3,5,9),4] <- 1
	A[c(4,6,9,10),5] <- 1
	A[c(1,5,7,10),6] <- 1
	A[c(1,2,6,8,10),7] <- 1
	A[c(2,3,7,9),8] <- 1
	A[c(3,4,5,8,10),9] <- 1
	A[c(5,6,7,9),10] <- 1

	g <- as.graph(A)
   g$plot.mode="COORDS"

	coord1 <- cbind(2*sin(2 * pi * ((0:(7 - 2))/(7-1))),
						 2*cos(2 * pi * ((0:(7 - 2))/(7-1))))
	coord2 <- rbind(c(0,1),c(1,0),c(0,-1),c(-1,0))
	g$coord=rbind(coord1,coord2)
	g
}

ladder <- function(rungs,open=F)
{
	 if(rungs<3) stop("must have at least 3 rungs")
    n <- 2*rungs
	 A <- matrix(0,nrow=n,ncol=n)
	 g1 <- chain(rungs,gmode="graph")
	 A[1:rungs,1:rungs] <- adjacency.matrix(g1)
	 A[(rungs+1):n,(rungs+1):n] <- adjacency.matrix(g1)
	 if(open){
	    start <- 2
		 end <- rungs-1
	 }
	 else{
	    start <- 1
		 end <- rungs
	 }
	 for(i in start:end){
	    A[i,rungs+i] <- 1
	    A[rungs+i,i] <- 1
	 }
	 g <- as.graph(A)
	 g$plot.mode <- "COORDS"
	 g$coord <- rbind(cbind(1:rungs,rep(0,rungs)),cbind(1:rungs,rep(1,rungs)))
	 g
}

heawood <- function()
{
   g <- Cn(14,gmode="graph")
	A <- adjacency.matrix(g)
	A[1,10] <- 1
	A[2,7] <- 1
	A[3,12] <- 1
	A[4,9] <- 1
	A[5,14] <- 1
	A[6,11] <- 1
	A[8,13] <- 1
	A[10,1] <- 1
	A[7,2] <- 1
	A[12,3] <- 1
	A[9,4] <- 1
	A[14,5] <- 1
	A[11,6] <- 1
	A[13,8] <- 1
	g$A <- A
	g
}

coxeter <- function()
{
   n <- 28
	A <- matrix(0,nrow=n,ncol=n)
	A[1,c(2,7,9)] <- 1
	A[2,c(1,3,10)] <- 1
	A[3,c(2,4,11)] <- 1
	A[4,c(3,5,12)] <- 1
	A[5,c(4,6,13)] <- 1
	A[6,c(5,7,14)] <- 1
	A[7,c(1,6,8)] <- 1
	A[8,c(15,28)] <- 1
	A[9,c(16,22)] <- 1
	A[10,c(17,23)] <- 1
	A[11,c(18,24)] <- 1
	A[12,c(19,25)] <- 1
	A[13,c(20,26)] <- 1
	A[14,c(21,27)] <- 1
	A[15,c(17,20)] <- 1
	A[16,c(18,21)] <- 1
	A[17,19] <- 1
	A[18,20] <- 1
	A[19,21] <- 1
	A[22,c(25,26)] <- 1
	A[23,c(26,27)] <- 1
	A[24,c(27,28)] <- 1
	A[25,c(22,28)] <- 1
	g <- symmetrize.graph(as.graph(A))
	g$plot.mode <- "COORDS"

	coord4 <- cbind(sin(2 * pi * ((0:(7 - 1))/(7))),
						  cos(2 * pi * ((0:(7 - 1))/(7))))
	coord3 <- cbind(2*sin(2 * pi * ((0:(7 - 1))/(7))-.6),
						 2*cos(2 * pi * ((0:(7 - 1))/(7))-.6))
	coord2 <- cbind(3*sin(2 * pi * ((0:(7 - 1))/(7))-.4),
						 3*cos(2 * pi * ((0:(7 - 1))/(7))-.4))
	coord1 <- cbind(4*sin(2 * pi * ((0:(7 - 1))/(7))),
						 4*cos(2 * pi * ((0:(7 - 1))/(7))))

   g$coord <- rbind(coord1,coord2,coord3,coord4)
	g
}

prism <- function(n)
{
	if(n<3) stop("n must be at least 3")
	g <- Cn(n,gmode="graph")
	p <- graph.adjoin(g,g,cbind(1:n,1:n))
	p$plot.mode <- "mcircles"
	p$m <- 2
	p
}

sun <- function(n)
{
	if(n<3) stop("n must be at least 3")
	 g <- Cn(n)
	 h <- empty.graph(n)
	 G <- graph.adjoin(h,g,cbind(1:n,1:n))
	 G$plot.mode <- "mcircles"
	 G$m <- 2
	 G
}

thick.sun <- function(n)
{
	if(n<3) stop("n must be at least 3")
	 g <- Cn(n)
	 h <- empty.graph(n)
	 G <- graph.adjoin(h,g,rbind(cbind(1:n,1:n),cbind(1:n,c(n,1:(n-1)))))
	 G$plot.mode <- "mcircles"
	 G$rot <- T
	 G$m <- 2
	 G
}

star.tail <- function(n,m)
{
   g <- star(n)
	h <- chain(m,gmode="graph")
   G <- graph.adjoin(g,h,c(floor((n-1)/2),1))
	G$plot.mode <- "COORDS"
	G$coord <- rbind(g$coord,cbind(1:m+1,rep(g$coord[floor((n-1)/2),2],m)))
	G
}

square.lattice <- function(n,m)
{
   g <- chain(n,gmode="graph")
	coord <- cbind(1:n,rep(1,n))
	k <- 0
	for(i in 2:m){
	   g <- graph.adjoin(g,chain(n,gmode="graph"),cbind(1:n+k,1:n))
		coord <- rbind(coord,cbind(1:n,rep(i,n)))
		k <- n+k
	}
	g$plot.mode <- "COORDS"
	g$coord <- coord
	g
}

star.of.david <- function()
{
   g <- thick.star(6)
	g$plot.mode <- "COORDS"
	g$coord <- rbind(c(2,4),c(4,3),c(4,1),c(2,0),c(0,1),c(0,3),
	                 c(2.666,3),c(3.333,2),c(2.666,1),c(1.333,1),c(0.666,2),
						  c(1.333,3))
	g
}

dude <- function(head=16,body=14,arm=7,leg=7)
{
	n <- head
	g1 <- Cn(n)
   x1 <- cbind(.5*sin(pi+2 * pi * ((0:(n - 1)/(n)))),
              .5*cos(pi+2 * pi * ((0:(n - 1)/(n)))))
	n <- body
   g2 <- chain(n)
	x2 <- cbind(rep(0,n),seq(-1.5,-7.5,length=n))
	g <- graph.adjoin(g1,g2,c(1,1))
	x <- rbind(x1,x2)
	n <- arm
	q <- seq(-1,0,length=n+1)[1:n]
	x3 <- cbind(q,-2*(q+1)-2)
	g3 <- chain(n)
	g <- graph.adjoin(g,g3,c(head+body/2,n))
	x <- rbind(x,x3)
	x4 <- cbind(-q,-2*(q+1)-2)
	g <- graph.adjoin(g,g3,c(head+body/2,n))
	x <- rbind(x,x4)
	n <- leg
	q <- seq(-1,0,length=n+1)[1:n]
	x5 <- cbind(q,2*(q)-7.5)
	g4 <- chain(n)
	g <- graph.adjoin(g,g4,c(head+body,n))
	x <- rbind(x,x5)
	g <- graph.adjoin(g,g4,c(head+body,n))
	x6 <- cbind(-q,2*(q)-7.5)
	x <- rbind(x,x6)
	g$plot.mode <- "COORDS"
	g$coord <- x
	g$xlim <- range(as.vector(x))+4
	g$ylim <- range(as.vector(x))
	g
}

J <- function(v,k,i) 
{
	if(v<k || k<i || i<0) stop("must have v>=k>=i>=0")
   require(gregmisc)
	R <- combinations(v,k)
	n <- nrow(R)
	A <- matrix(0,nrow=n,ncol=n)
	for(j in 1:(n-1)){
	   for(l in (j+1):n){
		   if(length(intersect(R[j,],R[l,]))==i) A[j,l] <- 1
		}
	}
	as.graph(A)
}

johnson <- function(v,k)
{
   if(v<2*k) stop("v must be at least 2*k")
	J(v,k,k-1)
}

knesser <- function(v,k)
{
   if(v<2*k) stop("v must be at least 2*k")
   J(v,k,0)
}

circulant <- function(n,C,mode=NULL)
{
	if(is.null(mode))
		if(all(sort(C)==sort(-C))) mode <- "graph"
	else mode <- "digraph"
   v <- 0:(n-1)
	C <- C %% n
	A <- matrix(0,nrow=n,ncol=n)
	for(i in 1:n)
	   for(j in 1:n){
		   if((i != j) && (((j-i)%%n) %in% C)) 
			   A[i,j] <- 1
		}
	if(mode=="digraph")
		g <- as.digraph(A)
	else
		g <- as.graph(A)
	g
}

deterministic.scale.free <- function(g,iterations=10,compute.coords=F)
{
   if(missing(g)){
      g <- empty.graph(3)
	  g <- add.edge(g,c(1,1),2:3)
   }
   for(i in 1:iterations){
	  if(compute.coords){
		 x <- graph.coords(g)
		 xr <- range(x[,1])
		 xx <- 2.0*diff(xr)
		 yr <- range(x[,2])
		 yy <- 1.2*diff(yr)
		 y <- cbind(x[,1]-xx,x[,2]-yy)
		 z <- cbind(x[,1]+xx,x[,2]-yy)
	  }
      h <- graph.adjoin(g,g,cbind(rep(1,graph.order(g)),1:graph.order(g)))
      g <- graph.adjoin(h,g,cbind(rep(1,graph.order(h)),1:graph.order(g)))
	  if(compute.coords){
		 g$coord <- rbind(x,y,z)
		 g$plot.mode <- "COORDS"
	  }
   }
   g
}

everett <- function()
{
   g <- graph.adjoin(Kn(4),trivial.graph(),rbind(c(2,1),c(3,1)))
	g <- delete.edge(g,2,3)
	g <- graph.adjoin(g,trivial.graph(),c(5,1))
	g <- graph.adjoin(g,Kn(4),rbind(c(6,1),c(6,4)))
	g <- delete.edge(g,7,10)
	coord <- rbind(
	   c(-2,1),
	   c(1,1),
	   c(1,-1),
	   c(-2,-1),
	   c(2,0),
	   c(3,0),
	   c(4,1),
	   c(7,1),
	   c(7,-1),
	   c(4,-1))
		g$plot.mode <- "COORDS"
		g$coord <- coord
		g
}

independent.set <- function(g)
{
	if(is.digraph(g)) stop("not implemented for digraphs")
   G <- graph.complement(g)
	clique(G)
}

independence.number <- function(g)
{
   clique.number(graph.complement(g))
}
intersection.graph <- function(R, m=NULL, withprob=F,mode="graph",tiny=F,threshold=1)
{
	n <- length(R)
	A <- matrix(0,nrow=n,ncol=n)
	if(!withprob) mode <- "graph"
	else P <- matrix(0,nrow=n,ncol=n)
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			if(withprob){
				if(is.null(m)){
					P[i,j] <- length(intersect(R[[i]],R[[j]]))/m
					P[j,i] <- P[i,j]
				}
				else {
					P[i,j] <- length(intersect(R[[i]],R[[j]]))/max(R[[i]],R[[j]])
					P[j,i] <- P[i,j]
				}
				q <- rbinom(2,1,P[i,j])
				A[i,j] <- q[1]
				if(mode=="graph")
				   A[j,i] <- A[i,j]
			   else
				   A[j,i] <- q[1]
			}
			else {
				A[i,j] <- length(intersect(R[[i]],R[[j]]))>=threshold
				A[j,i] <- A[i,j]
			}
		}
	}
	if(tiny){
		if(mode=="graph")
			g <- as.graph(A)
		else 
			g <- as.digraph(A)
	}
	else{
		if(mode=="graph")
			g <- as.graph(A)
		else 
			g <- as.digraph(A)
		g$R <- R
		if(withprob)
			g$p <- p
		g$threshold <- threshold
		class(g) <- c("IG","graph")
	}
   g
}

sphere.intersection.graph <- function(x=NULL,R,d=NULL, p=2,tiny=F)
{
	if(is.null(d)){
	  if(is.null(x)) stop("one of x or d must be given")
	  d <- pdist(x,p=p)
	}
	if(missing(R)) stop("radii must be given")
	n <- nrow(d)
	if(length(R)< n) R <- rep(R,n)[1:n]
	A <- matrix(0,nrow=n,ncol=n)
	for(i in 1:n){
		A[i,] <- d[i,]<(R[i]+R)
	}

	if(tiny){
		g <- as.graph(A)
	}
	else{
		if(is.null(x)){
			g <- as.graph(A)
			g$R <- R
			g$d <- d
			g$p <- p
		}
		else{
			g <- as.graph(A)
			g$R <- R
			g$d <- d
			g$p <- p
			g$coord <- x
			g$plot.mode <- "COORDS"
		}
		class(g) <- c("SphereIG","graph")
	}
   g
}

coin.graph <- function(x=NULL,r=1,d=NULL,p=2,tiny=F)
{
	if(is.null(x)){
	   if(is.null(d)) stop("one of x or d must be given")
		n <- nrow(d)
	}
	else n <- nrow(x)
   R <- rep(r,n)
	if(is.null(x))
		sphere.intersection.graph(R=R,d=d,p=p,tiny=tiny)
   else
		sphere.intersection.graph(R=R,x=x,d=d,p=p,tiny=tiny)
}

common.IG <- function(g)
{
	if(is.graph(g))
		R <- g$R
	else R <- g
	if(is.null(R) | !is.list(R)) stop("g must contain sets (R) or be a list of sets")
	a <- NULL
	A <- g$A
	a <- R[[1]]
	B <- R
	for(j in 1:length(R)) A[1,j] <- length(intersect(R[[1]],R[[j]]))
	for(i in 2:length(R)){
		a <- intersect(A,R[[i]])
		for(j in 1:length(R)){
		   A[i,j] <- length(intersect(R[[i]],R[[j]]))
		   B[[i]] <- setdiff(R[[i]],unlist(R[-i]))
		}
	}
	list(all=a,pairs=A,unique=B)
}

is.IG <- function(g,checkR=T)
{
   a <- is.graph(g) & inherits(g,"IG")
	if(checkR){
		a <- a & !is.null(g$R) & is.list(g$R)
	}
	a
}

test.intersect <- function(g,FUN=clique.number,n=100,threshold=g$threshold,side="upper",
                           nullstr="null",verbose=T,sets)
{
	if(!is.IG(g,checkR=T)) stop("Must be an intersection graph")
	cn <- FUN(g)
   if(is.null(g$R)) stop("Intersection graph must contain sets")
	if(graph.size(g)==0) return(rep(1,n))
	if(g$threshold != threshold){
	   g <- intersection.graph(R=g$R,threshold=threshold)
	}
	Q <- lapply(g$R,as.character)
	R <- g$R[unlist(lapply(Q,function(l)pmatch(nullstr,l,nomatch=F)!=1))]
	m <- length(R)
	if(missing(sets))
		sets <- unique(unlist(R))
	numbers <- unlist(lapply(R,length))
	clqs <- rep(0,n)
	for(i in 1:n){
	   r <- list(0)
		for(j in 1:m){
		   r[[j]] <- sample(sets,numbers[j])
		}
		G <- intersection.graph(R=r)
		clqs[i] <- FUN(G)
	}
	if(verbose){
		cat("\nOriginal Value:",cn,"\n")
		cat("Number Vertices:",m,"\n")
		cat(n,"Monte Carlo runs\n")
		if(side=="upper")
			cat("Upper Pvalue:",sum(clqs>=cn)/n,"\n\n")
		else if(side=="lower")
			cat("Lower Pvalue:",sum(clqs<=cn)/n,"\n\n")
		else
			cat("Two-sided Pvalue:",sum(abs(clqs-mean(clqs))>abs(mean(clqs)-cn))/n,"\n\n")
	}
	invisible(clqs)
}


interval.graph <- function(intervals)
{
   n <- nrow(intervals)
	m <- apply(intervals,1,mean)
	R <- matrix(intervals[,2]-m,ncol=1)
	g <- sphere.intersection.graph(x=m,R=R)
	class(g) <- c(class(g),"Interval")
	g
}

juggle1 <- function(data,classes,sampled=T,
                   num=100,method="greedy",sample.proportion=0.1)
{
	if(sample.proportion<1.0) sampled<-T
	replace <- sample.proportion>=1
	cls <- 1:length(unique(classes))
	if(is.null(data)) stop("data must be given")
	if(is.matrix(data)){
		d <- ncol(data)
	}
	else d <- 1
	dxx <- pdist(data,p=2,d=d)
	n <- nrow(dxx)
	S <- list(0)
	R <- list(0)
	for(i in 1:num){
		S[[i]] <- list(0)
		R[[i]] <- list(0)
		xindex <- list(0)
		for(j in cls){
			xindex[[j]] <- which(classes==j)
			if(sampled){
				nx <- sample.proportion*length(xindex[[j]])
				xindex[[j]] <- sort(unique(sample(xindex[[j]],nx,replace=replace)))
			}
		}
		for(j in cls){
			DXX <- dxx[xindex[[j]],xindex[[j]]]
			ind <- setdiff(unlist(xindex),xindex[[j]])
			DYX <- dxx[ind,xindex[[j]]]
			rx <- apply(DYX,2,min)
			Gx <- as.digraph(matrix(as.integer(DXX<rx),nrow=nrow(DXX)))
			sx <- dominate.random.sample(Gx,method)
			S[[i]][[j]] <- xindex[[j]][sx]
			R[[i]][[j]] <- rx[sx]
		}
	}
	list(S=S,R=R,dimension=d)
}

juggle1.classify <- function(data,J,tdata,indices)
{
	if(missing(indices)) indices <- 1:length(J$S)
	n <- length(indices)
	if(any(indices>length(J$S))) stop("invalid indices")
	if(J$dimension>1)
		N <- nrow(data)
	else N <- length(data)
	nc <- length(J$S[[1]])
	out <- matrix(0,nrow=N,ncol=nc)
	inds <- sort(unique(unlist(J$S)))
	if(J$dimension>1)
		d <- pdist(data,tdata[inds,],d=J$dimension)
	else
		d <- pdist(data,tdata[inds],d=1)
	foo <- matrix(0,nrow=nc,ncol=N)
	for(i in indices){
		for(j in 1:nc){
			if(length(J$S[[i]][[j]])>1)
				if(N==1)
					foo[j,] <- min(d[,match(J$S[[i]][[j]],inds)]/J$R[[i]][[j]])
				else
					foo[j,] <- apply(t(t(d[,match(J$S[[i]][[j]],inds)])/J$R[[i]][[j]]),1,min)
			else
				foo[j,] <- d[,match(J$S[[i]][[j]],inds)]/J$R[[i]][[j]]
		}
		a <- apply(foo,2,smallestplace)
		for(j in 1:nc){
			out[,j] <- out[,j]+1*(a==j)
		}
	}
	out/n
}

juggle2 <- function(data,classes,sampled=T,
                   num=100,method="greedy",sample.proportion=0.1,
						 k=2)
{
	if(sample.proportion<1.0) sampled<-T
	replace <- sample.proportion>=1
	cls <- 1:length(unique(classes))
	if(is.null(data)) stop("data must be given")
	if(is.matrix(data)){
		d <- ncol(data)
	}
	else d <- 1
	if(k>=d) stop("use juggle if no dimension sampling is to be performed")
	n <- nrow(data)
	S <- list(0)
	R <- list(0)
	vars <- list(0)
	for(i in 1:num){
		S[[i]] <- list(0)
		R[[i]] <- list(0)
		xindex <- list(0)
		for(j in cls){
			xindex[[j]] <- which(classes==j)
			if(sampled){
				nx <- sample.proportion*length(xindex[[j]])
				xindex[[j]] <- sort(unique(sample(xindex[[j]],nx,replace=replace)))
			}
		}
		if(k>0) K<-k
		else{
		   K <- rbinom(1,d,.5)
		}
		vars[[i]] <- sample(1:d,K)
		dxx <- pdist(data[,vars[[i]]],p=2,d=K)
		for(j in cls){
			DXX <- dxx[xindex[[j]],xindex[[j]]]
			ind <- setdiff(unlist(xindex),xindex[[j]])
			DYX <- dxx[ind,xindex[[j]]]
			rx <- apply(DYX,2,min)
			Gx <- as.digraph(matrix(as.integer(DXX<rx),nrow=nrow(DXX)))
			sx <- dominate.random.sample(Gx,method)
			S[[i]][[j]] <- xindex[[j]][sx]
			R[[i]][[j]] <- rx[sx]
		}
	}
	list(S=S,R=R,dimension=d,vars=vars)
}

juggle2.classify <- function(data,J,tdata,indices)
{
	if(missing(indices)) indices <- 1:length(J$S)
	n <- length(indices)
	if(any(indices>length(J$S))) stop("invalid indices")
	if(J$dimension>1)
		N <- nrow(data)
	else N <- length(data)
	nc <- length(J$S[[1]])
	out <- matrix(0,nrow=N,ncol=nc)
	inds <- sort(unique(unlist(J$S)))
	foo <- matrix(0,nrow=nc,ncol=N)
	for(i in indices){
		d <- pdist(data[,J$vars[[i]]],tdata[inds,J$vars[[i]]],d=length(J$vars[[i]]))
		for(j in 1:nc){
			if(length(J$S[[i]][[j]])>1)
				if(N==1)
					foo[j,] <- min(d[,match(J$S[[i]][[j]],inds)]/J$R[[i]][[j]])
				else
					foo[j,] <- apply(t(t(d[,match(J$S[[i]][[j]],inds)])/J$R[[i]][[j]]),1,min)
			else
				foo[j,] <- d[,match(J$S[[i]][[j]],inds)]/J$R[[i]][[j]]
		}
		a <- apply(foo,2,smallestplace)
		for(j in 1:nc){
			out[,j] <- out[,j]+1*(a==j)
		}
	}
	out/n
}


juggle <- function(data,classes,sampled=T,sample.dim=F,
                   num=100,method="greedy",sample.proportion=0.1,
						 k=2)
{
   if(sample.dim) a <- juggle2(data,classes,sampled,num,
	                            method,sample.proportion,k)
   else a <- juggle1(data,classes,sampled,num,method,sample.proportion) 
	a
}

juggle.classify <- function(data,J,tdata,indices)
{
   if(is.null(J$vars)) a <- juggle1.classify(data,J,tdata,indices)
   else a <- juggle2.classify(data,J,tdata,indices)
	a
}
L <- function(g){
   diag(degrees(g,cmode="indegree"))-adjacency.matrix(g)
}  

laplacian <- function(g){
	d <- degrees(g,cmode="indegree")
   Tt <- 1/sqrt(d)
	Tt[d==0] <- 0
	Tt <- diag(Tt)
	Tt %*% L(g) %*% Tt
}

seidel <- function(g){
   A <- -adjacency.matrix(g) + 1*!adjacency.matrix(g)
	diag(A) <- 0
	A
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

mst <- function(x=NULL,distmat=NULL, p=2,tiny=F)
{
	if(is.graph(x)){
		g <- x
	   if(!is.connected(g)) stop("the graph must be connected")
		x <- graph.coords(g)
		distmat <- geodesic.dist(g)
	}
	if(is.null(distmat)){
		if(is.null(x)) stop("One of x or distmat must be given")
		if(is.matrix(x)){
			distmat <- pdist(x,p=p)
		}
		else if(is.graph(x,ignore.symmetry=F)){
		   distmat <- geodesic.dist(x)
		}
		else if(is.digraph(x)){
		   warning("converting x to a graph")
			distmat <- geodesic.dist(as.graph(x,symmetric=T))
		}
	}
   n <- nrow(distmat);
   retval <- .C("mst",
				graph  = as.double(distmat),
				mst    = integer(2*n),
				n      = as.integer(n),
				nedges = as.integer(0),
				sum    = as.double(0))
   mst <- matrix(retval$mst,byrow=T,nrow=nrow(distmat))
   nedges <- retval$nedges
   mst <- mst[1:nedges,]+1
   sum <- retval$sum
	A <- matrix(0,ncol=ncol(distmat),nrow=nrow(distmat))
	for(i in 1:nedges) {
		A[mst[i,1],mst[i,2]] <- 1
		A[mst[i,2],mst[i,1]] <- 1
	}
	if(tiny){
		out <- as.graph(A)
	}
	else{
		out <- as.graph(A)
		out$edges <- mst
		out$nedges <- nedges
		out$sum <- sum
		out$coord <- x
		out$plot.mode <- "COORDS"
		class(out) <- c("MST","graph")
	}
	out
}

number.spanning.trees <- function(g)
{
   det(L(g)[-1,-1])
}
nng <- function(x=NULL,dx=NULL,k=1,mutual=F,p=2,tiny=F,mode="digraph")
{
   if(is.null(dx)) {
	  if(is.null(x)) stop("one of x or dx must be given")
	  dx <- pdist(x,p=p)
   }
   n <- nrow(dx)
   A <- matrix(0,nrow=n,ncol=n)
   for(i in 1:n){
	  d <- sort(dx[i,])
      A[i,dx[i,]<=d[k+1]] <- 1
   }
   diag(A) <- 0
   if(mutual){
      for(i in 1:n){
	     A[i,] <- A[i,] & A[,i]
		 A[,i] <- A[i,]
	  }
   }
	if(tiny){
		if(mutual)
		   out <- as.graph(A)
		else
		   out <- as.digraph(A)
	}
	else{
		if(mutual)
		   out <- as.graph(A)
		else
		   out <- as.digraph(A)
		out$k <- k
		out$mutual <- mutual
		out$p <- p
		out$dx <- dx
		out$coord <- x
		out$plot.mode <- "COORDS"
		class(out) <- c("NNG","graph")
	}
	m <- charmatch(mode,c("graph","digraph"))
	if(m==1 && !mutual) out <- as.graph(symmetrize.graph(out))
   out
}

nodal.domains <- function(G,k=NULL,type="weak",epsilon=10^(-10))
{
	x <- pmatch(type,c("strong","weak"))
	if(is.na(x)) stop("Could not disambiguate type: must be \"weak\" or \"strong\"")
	n <- graph.order(G)
	if(!is.null(k)) {
		if(k>n) stop("k must be less than the |V|")
		if(k<=0) stop("k must be positive")
	}
	comp <- graph.components(G)
	f <- fiedler(G)
	if(is.null(k)) k <- f$index
	lambda <- f$values[k]
	mult <- sum(abs(f$values-lambda)<epsilon)
	nodal.domains <- list(0)
	nodal.signs <- list(0)
	ks <- (1:n)[abs(f$values-lambda)<epsilon]
	for(l in 1:mult){
		v <- f$vectors[,ks[l]]
		y <- max(comp)
		d <- comp
		v[abs(v)<epsilon] <- 0

		d[v==0] <- 0
		d[v>0] <- d[v>0]+y

		signs <- NULL

		domains <- list(0)
		for(i in 1:max(d)){
			if(!any(d==i)){
				domains[[i]] <- 0
				signs[i] <- 0
			}
			else{
				domains[[i]] <- (1:n)[d==i]
				signs[i] <- sign(v[d==i][1])
			}
		}

		if(x==2){
			for(i in 1:n){
				if(d[i]==0){
					for(j in 1:length(domains)){
						domains[[j]] <- setdiff(c(domains[[j]],i),0)
					}
				}
			}
		}
		nodal.domains[[l]] <- domains
		nodal.signs[[l]] <- signs
	}
	list(domains=nodal.domains,
	     signs=nodal.signs,
	     v=f$vectors[,k:(k+mult-1)],
		  k=ks,
		  value=lambda,
		  multiplicity=mult)
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

guess.dim <- function(x,y)
{
   if(is.vector(x)){
      if(is.null(y)){
	     d <- 1
	  }
	  else{
	     if(is.matrix(y)){
		    d <- ncol(y)
		 }
		 else{
		    if(length(x)==length(y)){
			   d <- length(x)
			}
			else{
			   d <- 1
			}
		 }
	  }
   }
   else{
      d <- ncol(x)
   }
   d
}

pdist <- function(x,y=NULL,d,p=2,w=NULL)
{
   if(missing(d)) d <- guess.dim(x,y)
   x <- matrix(x,ncol=d)
   nx <- nrow(x)
   if(!is.null(y)){
	  y <- matrix(y,ncol=d)
	  return(pdistxy(x,y,d,p,w))
   }
   out <- rep(0,nx^2)
   if(is.infinite(p)){
	  retval <- .C("pdistinf",
				   x = as.double(t(x)),
				   n = as.integer(nx),
				   d = as.integer(d),
				   w = as.double(w),
				   out=as.double(out))
   }
   else {
	  retval <- .C("pdist",
				   x = as.double(t(x)),
				   n = as.integer(nx),
				   d = as.integer(d),
				   p = as.double(p),
				   w = as.double(w),
				   out=as.double(out))
   }
   return(matrix(retval$out,nrow=nrow(x)))
}

pdistxy <- function(x,y,d,p=2,w=NULL)
{
   if(missing(d)){
	  d <- guess.dim(x,y)
   }
   x <- matrix(x,ncol=d)
   nx <- nrow(x)
   y <- matrix(y,ncol=d)
   ny <- nrow(y)

   out <- rep(0,nx*ny)
   if(is.infinite(p)){
	  retval <- .C("pdistxyinf",
				   x = as.double(t(x)),
				   y = as.double(t(y)),
				   nx = as.integer(nx),
				   ny = as.integer(ny),
				   d = as.integer(d),
				   w = as.double(w),
				   dist=as.double(out))
   }
   else {
	  retval <- .C("pdistxy",
				   x = as.double(t(x)),
				   y = as.double(t(y)),
				   nx = as.integer(nx),
				   ny = as.integer(ny),
				   d = as.integer(d),
				   p = as.double(p),
				   w = as.double(w),
				   dist=as.double(out))
   }
   matrix(retval$dist,byrow=T,nrow=nx)
}

pinwheel <- function(classes, R=length(unique(classes)),r=1,
				center=NULL)
{
   u <- unique(classes)
   coord <- matrix(0,nrow=length(classes),ncol=2)
   nc <- length(u)
	if(length(r)==1) r <- rep(r,nc)
	if(is.null(center)) {
		center <- -10^10
		theta <- seq(0,2*pi,length=nc+1)
	}
	else{
		if(any(classes==center))
			theta <- seq(0,2*pi,length=nc)
		else
			theta <- seq(0,2*pi,length=nc+1)
	}
	j <- 1
   for(i in 1:nc){
      m <- sum(classes==u[i])
      phi <- seq(0,2*pi,length=m+1)[1:m]
      cr <- R+r[i]
      x1 <- c(cr*cos(theta[j]),cr*sin(theta[j]))
		if(u[i]==center){
			coord[classes==u[i],] <- cbind(r[i]*cos(phi),r[i]*sin(phi))
		}
		else {
			coord[classes==u[i],] <- cbind(r[i]*cos(phi)+x1[1],r[i]*sin(phi)+x1[2])
			j <- j+1
		}
   }
	coord
}
midarrows <- function(x0, y0, x1, y1, length = 0.25, angle = 30, code = 2,
            col = par("fg"), lty = NULL, lwd = par("lwd"), xpd = FALSE)
{
   if(code==1){
      arrows(x1, y1, (x1+x0)/2, (y1+y0)/2, length, angle, code = 2,
	              col, lty, lwd, xpd)
      segments((x1+x0)/2, (y1+y0)/2, x0, y0, 
	              col, lty, lwd, xpd=xpd)
   }
   else if(code==2){
      arrows(x0, y0, (x1+x0)/2, (y1+y0)/2, length, angle, code = 2,
	              col, lty, lwd, xpd)
      segments((x1+x0)/2, (y1+y0)/2, x1, y1, 
	              col, lty, lwd, xpd=xpd)
   }
   else if(code==3){
      arrows(x0, y0, (x1+x0)/2, (y1+y0)/2, length, angle, code = 2,
	              col, lty, lwd, xpd)
      arrows(x1,y1,(x1+x0)/2, (y1+y0)/2, length, angle, code = 2,
	              col, lty, lwd, xpd)
   }
}

endarrows <- function(x0, y0, x1, y1, length = 0.05, code = 2,
				lty=1,lwd=par("lwd"),
            hcol=2,hlwd=1,col = par("fg"), ...)
{
   segments(x0,y0,x1,y1,col=col,lty=lty,lwd=lwd,...)
   if(code==1){
	  segments(x0,y0,x0+length*(x1-x0),y0+length*(y1-y0),col=hcol,lwd=hlwd,lty=lty,...)
   }
   else if(code==2){
	  segments(x1-(length)*(x1-x0),y1-(length)*(y1-y0),x1,y1,col=hcol,lwd=hlwd,lty=lty,...)
   }
   else if(code==3){
	  segments(x0,y0,x0+length*(x1-x0),y0+length*(y1-y0),col=hcol,lwd=hlwd,lty=lty,...)
	  segments(x1-(length)*(x1-x0),y1-(length)*(y1-y0),x1,y1,col=hcol,lwd=hlwd,lty=lty,...)
   }
}

plot.graph <- function (x, gmode = NULL, 
    label,label.coord=NULL,
    coord = NULL, jitter = FALSE, usearrows = TRUE, 
    mode = "circle", displayisolates = TRUE, pad = 0, vertex.pch = 19, 
    label.cex = 1, vertex.cex = 1, label.col = 1, edge.col = 1, 
    vertex.col = 1, arrowhead.length = 0.2, edge.type = 1, edge.lwd = 0, 
	 xlab=NULL,ylab=NULL,
	 newplot=T, arrowtype="arrow",head.lty=edge.type,head.lwd=edge.lwd,
	 head.col=edge.col+1,xlim,ylim,plotit=T,m=2,rot=F,
	 clique.cex=2,clique.label.cex=label.cex[1],clique.label.col=label.col[1],
	 keep.circle=F,clique.num=NULL,shadow=F,
	 cspring=1,crep=1,l=1,iterations=10,delta=0.01,start="random",
	 parts=NULL,force.mode=F,
	 bycomponent=F,component=0, 
	 R=3, r=1,center=NULL,
	 foffset=0,
	 ...)
{
	 g <- x
	 if(bycomponent){
	    g <- graph.component(g,component)
	 }
	 if(force.mode) g$plot.mode <- NULL
	 if(!force.mode & !is.null(g$plot.mode) & is.null(coord)){
		 mode <- g$plot.mode
		 if(!is.null(g$xlim)) xlim <- g$xlim
		 if(!is.null(g$ylim)) ylim <- g$ylim
		 if(mode=="bipartite") parts <- g$parts
		 if(mode=="pinwheel") parts <- g$parts
		 if(mode=="mpoly") parts <- g$parts
		 if(mode=="mflower") parts <- g$parts
		 else if(mode=="mcircles"){
			 m <- g$m
			 if(!is.null(g$rot)) rot <- g$rot
		 }
		 else if(mode=="COORDS"){
			 coord <- g$coord
			 if(is.null(coord))
			    mode <- "circle"
			 else
				 mode <- NULL
			 if(is.vector(coord)){
			    coord <- cbind(coord,runif(length(coord)))
			 }
		 }
	 }
	 if(!is.graph(g,ignore.symmetry=T) && !is.digraph(g)){
		 stop("object is not a graph")
	 }
	 if(missing(label)) label <- 1:graph.order(g)
	 ARROWS = c("arrow","midarrow","endarrow")
	 arrowtype <- pmatch(arrowtype,ARROWS)
	 if(arrowtype<1 || arrowtype>3){
	    warning("unknown arrowtype. using type \"arrow\"")
		 arrowtype<-1
	 }
	 d <- g$A
	 if(is.null(gmode)){
	    if(is.graph(g)) gmode <- "graph"
		 else gmode <- "digraph"
	 }
    if (gmode == "graph") {
        usearrows <- FALSE
    }
    n <- nrow(d)
    d[is.na(d)] <- 0
	 if(!is.null(mode)){
		 MODES <- c("princoord","eigen","seidel","laplacian","fiedler","mds","random","circle",
						"circrand","rmds","geodist","mst","dot","adj","seham","mcircles",
						"spring","tree","clique","bipartite","line","pinwheel","mpoly","mflower")
		 tmp <- charmatch(mode,MODES)
		 if(is.null(tmp)){
			 stop("invalid mode or mode not recognized")
		 }
		 else if(is.na(tmp)){
			 stop("invalid mode or mode not recognized")
		 }
		 else if(tmp==0){
			 stop("ambiguous mode: retry with more characters")
		 }
		 else if(tmp<1 || tmp>length(MODES)){
			 stop("invalid mode or mode not recognized")
		 }
		 mode <- MODES[tmp]
		 if(mode %in% 
			 c("princoord","eigen","seidel","laplacian","fiedler","adj","mds","sehan")){
				 if(is.null(xlab))
					 xlab <- expression(lambda[1])
				 if(is.null(ylab))
					 ylab <- expression(lambda[2])
			 }
		 else{
			 if(is.null(xlab))
				 xlab <- ""
			 if(is.null(ylab))
				 ylab <- ""
		 }
	 }
    if (!is.null(coord)) {
        x <- coord[, 1]
        y <- coord[, 2]
    }
	 else if (mode == "pinwheel"){
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
	    coord <- pinwheel(classes=parts,R=R,r=r,center=center)
            
		 x <- coord[,1]
		 y <- coord[,2]
	 }
	 else if (mode == "mpoly"){
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
	    coord <- mpoly(classes=parts)
		 x <- coord[,1]
		 y <- coord[,2]
	 }
	 else if (mode == "mflower"){
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
	    coord <- mflower(classes=parts)
		 x <- coord[,1]
		 y <- coord[,2]
	 }
    else if (mode == "princoord") {
        cd <- cor(rbind(d, t(d)), use = "pairwise.complete.obs")
        cd <- replace(cd, is.na(cd), 0)
        e <- eigen(cd, symmetric = TRUE)
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "eigen") {
        e <- eigen(d, symmetric = (gmode != "digraph"))
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "laplacian") {
        e <- eigen(laplacian(g), symmetric = (gmode != "digraph"))
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "fiedler") {
        e <- fiedler(g)
        x <- Re(e$vectors[, e$index+foffset])
        y <- Re(e$vectors[, e$index+1+foffset])
    }
    else if (mode == "seidel") {
        e <- eigen(seidel(g), symmetric = (gmode != "digraph"))
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "mds") {
        Dmat <- matrix(nrow = n, ncol = n)
        diag(Dmat) <- 0
        for (i in 1:n) for (j in 1:n) if (i > j) 
            Dmat[i, j] <- sqrt(sum(abs(d[i, ] - d[j, ])) + sum(abs(d[, 
                i] - d[, j])))
        Dmat[upper.tri(Dmat)] <- t(Dmat)[upper.tri(Dmat)]
        Imat <- matrix(nrow = n, ncol = n)
        Imat[, ] <- 0
        diag(Imat) <- 1
        Hmat <- Imat - (1/n) * rep(1, n) %*% t(rep(1, n))
        Amat <- -0.5 * Dmat^2
        Bmat <- Hmat %*% Amat %*% Hmat
        e <- eigen(Bmat)
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "random") {
        x <- runif(n, -1, 1)
        y <- runif(n, -1, 1)
    }
    else if (mode == "line") {
        x <- 1:n
        y <- rep(0,n)
    }
    else if (mode == "circle") {
        x <- sin(2 * pi * ((0:(n - 1))/n))
        y <- cos(2 * pi * ((0:(n - 1))/n))
    }
    else if (mode == "circrand") {
        tempd <- rnorm(n, 1, 0.25)
        tempa <- runif(n, 0, 2 * pi)
        x <- tempd * sin(tempa)
        y <- tempd * cos(tempa)
    }
    else if (mode == "rmds") {
        require(mva)
        tempmds <- cmdscale(dist(d))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "geodist") {
        require(mva)
        tempmds <- cmdscale(as.dist(geodesic.dist(g)))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "mst") {
		  q <- geodesic.dist(g)
        temp <- mst(distmat=q)
		  z <- dotlayout(temp)
        x <- z[, 1]
        y <- z[, 2]
    }
    else if (mode == "dot") {
		  z <- dotlayout(g)
        x <- z[, 1]
        y <- z[, 2]
    }
    else if (mode == "adj") {
        require(mva)
        tempmds <- cmdscale(as.dist(-d + max(d)))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "seham") {
        require(mva)
        require(sna)
        temp <- sedist(d)
        tempmds <- cmdscale(as.dist(temp))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "mcircles") {
		  if(m<2 || m>(n/2)) stop("invalid choice of m in mcircles")
		  r <- m:1
		  x <- NULL
		  y <- NULL
		  N <- n/m
		  for(i in 1:m){
			  x <- c(x,r[i]*sin(2 * pi * ((0:(N - 1))/N)-rot*(i %% 2)*pi/N))
			  y <- c(y,r[i]*cos(2 * pi * ((0:(N - 1))/N)-rot*(i %% 2)*pi/N))
		  }
		  N <- n-length(x)
		  if(N>0){
			  x <- c(x,r[m]*sin(2 * pi * ((0:(N - 1))/N)-rot*((m+1) %% 2)*pi/N))
			  y <- c(y,r[i]*cos(2 * pi * ((0:(N - 1))/N)-rot*((m+1) %% 2)*pi/N))
		  }
    }
	 else if (mode=="spring"){
		if(is.matrix(start)){
			if(nrow(start) != n || ncol(start) != 2) stop("start must be nx2")
			points <- start
		}
		else if(is.character(start)){
			points <- Recall(g, gmode = gmode, mode="random",plotit=F)
	 }
	 else stop("start must be either a mode name or a matrix")
		frep <- function(x,y,crep){
			dd <- pdist(x,y,d=2)
			crep/dd^3 * (x-y)
		}
		fspring <- function(x,y,cspring,l){
			dd <- pdist(x,y,d=2)
			cspring*log(dd/l)*(x-y)
		}
		force <- function(A,j,points,crep,cspring,l){
			f <- 0
			for(i in setdiff(1:nrow(A),j)){
				if(A[i,j]==1){
					f <- f+fspring(points[i,],points[j,],cspring,l)
				}
				else{
					f <- f+frep(points[i,],points[j,],crep)
				}
			}
			f
		}
		direction <- function(x,y){
			#acos(x %*% y/(sqrt(sum(x %*% x)*sum(y %*% y))))
			atan2(x ,y)
		}
		B <- d
		diag(B) <- 0
		oforces <- matrix(0,nrow=n,ncol=2)
		deltas <- rep(0,n)
		for(i in 1:iterations){
			forces <- matrix(0,nrow=n,ncol=2)
			for(j in 1:n){
				forces[j,] <- force(B,j,points,crep,cspring,l)
				if(i>1)
					dd <- direction(forces[j,],oforces[j,])
				else 
					dd <- 1
				if((dd>0) && (dd<.01)) {
					deltas[j] <- delta*2
				}
				else if((dd<0) && (dd> -0.01)) {
					deltas[j] <- delta/2
				}
				else {
					deltas[j] <- delta
				}
				oforces[j,] <- forces[j,]
			}
			for(j in 1:n){
				points[j,] <- points[j,]+delta*forces[j,]
			}
			#plot.graph(d,coord=points)
		}
		x <- points[,1]
		y <- points[,2]
	 }
	 else if(mode=="tree"){
	    v <- 1
		 x <- rep(0,n)
		 y <- rep(0,n)
		 S <- 1:n
		 i <- 1
		 while(length(S)>0){
			 q <- length(v)
			 k <- -q
			 for(j in 1:length(v)){
				 a <- neighborhood(g,v[j],open=T,as.indices=T)
				 b <- setdiff(a,v)
				 if(length(b)>0){
					 S <- setdiff(S,a)
					 z <- seq(k,k+1,length=length(b))
					 v <- c(v,a)
					 x[b] <- z
					 y[b] <- i
				 }
				 k <- k + length(b)+1
			 }
			 i <- i+1
		 }
	 }
	 else if(mode=="clique"){
		 if(is.complete(g)){
		    stop("graph is complete")
		 }
		 a <- clique.number(g)
		 if(!is.null(clique.num)){
			 if(a<clique.num) 
			    stop("clique.num must be less than the clique number of the graph")
		    a <- clique.num
		 }
		 if(a < n && a>2){
			 g <- collapse(g,vertices=clique(g,a))
			 if(!is.null(label)){
			    if(shadow){
			       slabel <- label[g$collapsed]
				 }
			    label <- c(paste(label[setdiff(1:n,g$collapsed)]),
				            eval(substitute(
								     expression(K[k]),list(k=length(g$collapsed)))))
			 }
			 d <- g$A
			 N <- nrow(g$A)
			 ov <- vertex.cex
			 vertex.cex <- rep(ov,N)
			 vertex.cex[N] <- clique.cex+ov
			 if(length(vertex.pch)>1 & length(vertex.pch)<n){
			    vertex.pch <- rep(vertex.pch,n/length(vertex.pch) +1)[1:n]
			 }
			 if(length(vertex.pch)==n){
			    vertex.pch <- vertex.pch[setdiff(1:n,g$collapsed)]
				 vertex.pch <- c(vertex.pch,13)
			 }
			 else{
				 vertex.pch <- rep(vertex.pch,N)
				 vertex.pch[N] <- 13
			 }
			 if(length(label.col)>1 & length(label.col)<n){
			    label.col <- rep(label.col,n/length(label.col) +1)[1:n]
			 }
			 if(length(label.col)==n){
			    label.col <- label.col[setdiff(1:n,g$collapsed)]
				 label.col <- c(label.col,clique.label.col)
			 }
			 else{
				 label.col <- rep(label.col,N)
				 label.col[N] <- clique.label.col
			 }
			 if(length(label.cex)>1 & length(label.cex)<n){
			    label.cex <- rep(label.cex,n/length(label.cex) +1)[1:n]
			 }
			 if(length(label.cex)==n){
			    label.cex <- label.cex[setdiff(1:n,g$collapsed)]
				 label.cex <- c(label.cex,clique.label.cex)
			 }
			 else{
				 label.cex <- rep(label.cex,N)
				 label.cex[N] <- clique.label.cex
			 }
			 if(keep.circle){
				  a <- sin(2 * pi * ((0:(n - 1))/n))
				  b <- cos(2 * pi * ((0:(n - 1))/n))
				  x <- c(a[setdiff(1:n,g$collapsed)],0)
				  y <- c(b[setdiff(1:n,g$collapsed)],0)
				  if(shadow){
				     z <- a[g$collapsed]
				     w <- b[g$collapsed]
				  }
			 }
			 else{
				 m <- N-1
				 x <- sin(2 * pi * ((0:(m - 1))/m))
				 y <- cos(2 * pi * ((0:(m - 1))/m))
				 x <- c(x,0)
				 y <- c(y,0)
			 }
			 if(!is.null(label.coord)){
				 xl <- c(label.coord[setdiff(1:n,g$collapsed),1],0)
				 yl <- c(label.coord[setdiff(1:n,g$collapsed),2],0)
				 if(shadow){
					 z <- c(label.coord[g$collapsed,1],0)
					 w <- c(label.coord[g$collapsed,2],0)
				 }
			 }
			 n <- N
		 }
		 else {
		    stop("Either g is complete, or contains no cliques")
		 }
	 }
	 else if(mode=="bipartite"){
		 if(is.null(parts)){
			 G <- symmetrize.graph(g)
			 a <- 1
			 b <- neighborhood(G,a,as.indices=T,open=T)
			 while(T){
			    a <- unique(neighborhood(G,b,as.indices=T,open=T))
			    b <- unique(neighborhood(G,a,as.indices=T,open=T))
				 if(length(c(a,b))>=n) break
			 }
			 if(length(intersect(a,b))>0) stop("g is not bipartite")
			 parts <- rep(0,n)
			 parts[a] <- 1
			 parts[b] <- 2
		 }
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
		 if(nrow(d) != ncol(d)){
			 a <- nrow(d)
			 b <- ncol(d)
			 n <- a+b
		    A <- matrix(0,nrow=n,ncol=n)
			 for(i in 1:a){
			    for(j in 1:b){
				    if(d[i,j]){
					    A[i,a+j] <- d[i,j]
					    A[a+j,i] <- d[i,j]
					 }
				 }
			 }
			 d <- A
			 if(!is.null(label)){
			    if(length(label)<n) label <- 1:n
			 }
			 usearrows <- F
			 x <- c(rep(0,a),rep(1,b))
			 y <- c((1:a)/a,(1:b)/b)
		 }
		 else{
			 if(length(parts) != n)
             stop("parts must be a vector of length the number of vertices")
			 uparts <- sort(unique(parts))
			 mp <- max(uparts)
			 if(any(uparts != 1:mp))
			    stop("parts must be a vector containing values from 1 to max(parts)")
		    A <- matrix(0,nrow=length(parts),ncol=length(parts))
			 x <- NULL
			 y <- NULL
			 a <- rep(0,mp)
			 for(i in 1:mp){
			    a[i] <- sum(parts==i)
				 x <- c(x,rep(i,a[i]))
				 y <- c(y,(1:a[i])/a[i])
			 }
			 b <- c(0,cumsum(a))
			 for(i in 1:mp){
			    for(j in 1:mp){
				    if(i != j){
						 A[(b[i]+1):(b[i+1]),(b[j]+1):b[j+1]] <- d[parts==i,parts==j]
					 }
				 }
			 }
			 d <- A
		 }
	 }
    if (jitter) {
		  x <- jitter(x,factor=as.integer(jitter))
		  y <- jitter(y,factor=as.integer(jitter))
    }
	 if(is.null(label.coord)){
		 xl <- x
		 yl <- y
	 }
	 else if(mode != "clique"){
		 xl <- label.coord[,1]
		 yl <- label.coord[,2]
	 }
    use <- displayisolates | !isolated.vertices(g)
	 if(missing(xlim)){
		 xlim <- c(min(x[use]) - pad, max(x[use]) + pad)
		 if(shadow==T && mode=="clique" && keep.circle==T){
		    xlim <- c(min(xlim[1],z-pad),max(xlim[2],z+pad))
		 }
	 }
	 if(missing(ylim)){
		 ylim = c(min(y[use]) - pad, max(y[use]) + pad) 
		 if(shadow==T && mode=="clique" && keep.circle==T){
		    ylim <- c(min(ylim[1],w-pad),max(ylim[2],w+pad))
		 }
	 }
	 if(plotit){
		 if(newplot){
			 if ((length(x) > 0) & (!all(use == FALSE)) & (any(vertex.cex>0))) 
				  plot(x[use], y[use], xlim = xlim, ylim = ylim,
						type = "p", pch = vertex.pch, xlab = xlab,
						ylab = ylab, col = vertex.col, cex = vertex.cex, 
						...)
			 else plot(0, 0, type = "n", pch = vertex.pch, xlab = xlab,
				  ylab = ylab, col = vertex.col, cex = vertex.cex, 
			  ...)
		 }
		 px0 <- vector()
		 py0 <- vector()
		 px1 <- vector()
		 py1 <- vector()
		 e.lwd <- vector()
		 e.col <- vector()
		 e.type <- vector()
		 for (i in c(1:n)[use]) for (j in c(1:n)[use]) if (i != j & d[i,j]>0){
			  px0 <- c(px0, as.real(x[i]))
			  py0 <- c(py0, as.real(y[i]))
			  px1 <- c(px1, as.real(x[j]))
			  py1 <- c(py1, as.real(y[j]))
			  if(is.matrix(edge.lwd)){
			     e.lwd <- c(e.lwd,edge.lwd[i,j])
			  }
			  else if (edge.lwd > 0) 
					e.lwd <- c(e.lwd, edge.lwd)
			  else e.lwd <- c(e.lwd, 1)
			  if(is.matrix(edge.col)){
			     e.col <- c(e.col,edge.col[i,j])
			  }
			  else if (edge.col > 0) 
					e.col <- c(e.col, edge.col)
			  else e.col <- c(e.col, 1)
			  if(is.matrix(edge.type)){
			     e.type <- c(e.type,edge.type[i,j])
			  }
			  else if (edge.type > 0) 
					e.type <- c(e.type, edge.type)
			  else e.type <- c(e.type, 1)
		 }
		 if (usearrows & (length(px0) > 0)) 
			  switch(arrowtype,
			  arrows(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), length = arrowhead.length, angle = 15, 
					col = e.col, lty = e.type, lwd = e.lwd),
			  midarrows(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), length = arrowhead.length, angle = 15, 
					col = e.col, lty = e.type, lwd = e.lwd),
			  endarrows(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), length = arrowhead.length, 
					col = e.col, hcol=head.col, hlwd=head.lwd,
					lty = e.type, lwd = e.lwd))
		 else if (length(px0) > 0) 
			  segments(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), col = e.col, lty = e.type, 
					lwd = e.lwd)
		 if ((length(label) > 0) & (!all(use == FALSE))) 
                      for(i in 1:length(xl)){
			  text(xl[i], yl[i], label[i], pos = 1, cex = label.cex, 
					col = label.col)
                      }
		 if ((length(x) > 0) & (!all(use == FALSE)) & (any(vertex.cex>0))) 
			  points(x[use], y[use], pch = vertex.pch,col = vertex.col, 
			  cex = vertex.cex) 
		 if(mode=="clique" && keep.circle==T && shadow==T){
			 text(z,w,slabel,pos=1,cex=label.cex[1],col=gray(0.8))
		 }
	 }
    invisible(cbind(x,y))
}
power.graph <- function(g,k=2,matmult=T)
{
	if(matmult){
	   A <- adjacency.matrix(g)
		diag(A) <- 1
		B <- A
		if(is.null(k)){
		  h <- list(0)
		  h[[1]] <- as.graph(B)
		  for(k in 2:max(as.vector(d))){
			  B <- ((B %*% A)>0)+0
			  h[[k]] <- as.graph(B)
		  }
		}
		else {
			for(i in 2:k){
			   B <- ((B %*% A)>0)+0
			}
			h <- as.graph(B)
		}
	}
	else{
		d <- geodesic.dist(g)
		if(is.null(k)){
		  h <- list(0)
			for(k in 1:max(as.vector(d))){
			  h[[k]] <- as.graph(d<=k)
		  }
		}
		else
		  h <- as.graph(geodesic.dist(g)<=k)
	}
   h
}

vector.product.graph <- function(vects,r=1,w=rep(1,ncol(vects)),tiny=F)
{
   n <- nrow(vects)
	A <- matrix(0,nrow=n,ncol=n)
	P <- (vects %*% (w*t(vects)))^r
	diag(P) <- 0
	if(any(P>1) || any(P<0)) stop("vectors must dot between 0 and 1")
	g <- rg(n,p=P)
	if(!tiny) {
		g$v <- vects
		g$p <- P
	}
	class(g) <- c(class(g),"Dot")
	g
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

qhull <- function(matrx,
				  joggle     = FALSE,
				  summary    = FALSE,
				  total      = FALSE,
				  verify	 = FALSE,
				  silent     = TRUE)
{
    if (!missing(matrx))
    {
        ### data matrix ###
        x <- eval(matrx)
		d <- dim(x)
		n <- d[1]
		dimension <- d[2]
        if (sum(abs(x[!is.na(x)]) == Inf) > 0)
        {
          cat("Sorry, qhull can't handle Infs\n")
          return()
        }
		if (sum(is.na(x)) > 0)
		{
          cat("Sorry, qhull can't handle NAs\n")
          return()
		}
        dfile <- tempfile("unix")
		ofile <- tempfile("unix")
		y <- c(dimension,n)
        write(y, file = dfile)
		write(t(x),file=dfile,append=T,ncol=dimension)

		command <- paste("cat",dfile,"| qhull",sep=" ")

		if(joggle)
		   command <- paste(command,"QJ",sep=" ")

		if(verify)
		   command <- paste(command,"Tv",sep=" ")

		if(summary)
		   command <- paste(command,"s",sep=" ")

		if(total)
		   command <- paste(command,"FA",sep=" ")

		command <- paste(command,"Fx",sep=" ")

		if(silent) command <- paste(command,"Pp",sep=" ")

        command <- paste(command, "TO",ofile,sep=" ")

        system(command, FALSE)

		q <- scan(ofile,quiet=TRUE)
		if(length(q) == 0){
		   return(NULL);
		}
		index <- q[2:length(q)]+1

		#clean up after yourself
		system(paste("/bin/rm",dfile,sep=" "))
		system(paste("/bin/rm",ofile,sep=" "))

		return(index)
    }
    else
      cat("Matrix argument required\n")
    NULL
}

qhull.hyperplanes <- function(matrx,
				  joggle     = FALSE,
				  summary    = FALSE,
				  total      = FALSE,
				  verify	 = FALSE,
				  silent     = TRUE)
{
    if (!missing(matrx))
    {
        ### data matrix ###
        x <- eval(matrx)
		d <- dim(x)
		n <- d[1]
		dimension <- d[2]
        if (sum(abs(x[!is.na(x)]) == Inf) > 0)
        {
          cat("Sorry, qhull can't handle Infs\n")
          return()
        }
		if (sum(is.na(x)) > 0)
		{
          cat("Sorry, qhull can't handle NAs\n")
          return()
		}
        dfile <- tempfile("unix")
		ofile <- tempfile("unix")
		y <- c(dimension,n)
        write(y, file = dfile)
		write(t(x),file=dfile,append=T,ncol=dimension)

		command <- paste("cat",dfile,"| qhull",sep=" ")

		if(joggle)
		   command <- paste(command,"QJ",sep=" ")

		if(verify)
		   command <- paste(command,"Tv",sep=" ")

		if(summary)
		   command <- paste(command,"s",sep=" ")

		if(total)
		   command <- paste(command,"FA",sep=" ")

		if(silent) 
		   command <- paste(command,"Pp",sep=" ")

		command <- paste(command,"n",sep=" ")

        command <- paste(command, "TO",ofile,sep=" ")

        system(command, FALSE)

		q <- scan(ofile,quiet=TRUE)
		if(length(q) == 0){
		   return(NULL);
		}
		index <- matrix(q[3:length(q)],byrow=T,ncol=(dimension+1))

		#clean up after yourself
		system(paste("/bin/rm",dfile,sep=" "))
		system(paste("/bin/rm",ofile,sep=" "))

		return(index)
    }
    else
      cat("Matrix argument required\n")
    NULL
}

qhull.interior <- function(y,hyperplanes,data)
{
   if(missing(hyperplanes)){
	  if(missing(data)){
	     cat("\nMust have either hyperplanes or data defined")
		 return(NULL)
	  }
      hyperplanes <- qhull.hyperplanes(data)
   }
   if(is.vector(y)){
      zz <- c(y,1)
	  ww <- hyperplanes %*% zz
	  sum(ww>1E-6) == 0
   }
   else {
	  zz <- cbind(y,rep(1,dim(y)[1]))
	  ww <- hyperplanes %*% t(zz)
	  aa <- apply(ww,2,">",1E-6)
	  apply(aa,2,sum) == 0
   }
}
rg <- function(n,p=0.5,mode="graph")
{
	 if(any(p<0) || any(p>1)) stop("p must be a probability")
	 MODES <- c("graph","digraph")
	 tmp <- charmatch(mode,MODES)
	 if(is.null(tmp)){
	    stop("invalid mode or mode not recognized")
	 }
	 else if(is.na(tmp)){
	    stop("invalid mode or mode not recognized")
	 }
	 else if(tmp==0){
	    stop("ambiguous mode: retry with more characters")
	 }
	 if(is.matrix(p)){
		 A <- matrix(0,nrow=n,ncol=n)
		 for(i in 1:n)
			 for(j in 1:n)
				 A[i,j] <- rbinom(1,1,p[i,j])
	 }
	 else if(length(p)>1){
		 A <- matrix(0,nrow=n,ncol=n)
		 for(i in 1:n){
			 A[i,] <- rbinom(n,1,p[i])
		 }
	 }
	 else {
		 if(mode=="digraph")
			 A <- matrix(rbinom(n*n,1,p),ncol=n,nrow=n)
		 else{
			 A <- matrix(0,nrow=n,ncol=n)
			 A[lower.tri(A)] <- rbinom(n*(n-1)/2,1,p)
		 }
	 }
	 diag(A) <- 0
	 if(mode=="graph"){
		 A[upper.tri(A)] <- t(A)[upper.tri(A)]
	 }
	 if(mode=="graph")
		 g <- as.graph(A)
	 else
		 g <- as.digraph(A)
	 g
}

rgdegree <- function(n,degree=2,mode="digraph",dmode="indegree", scale=T)
{
	 MODES <- c("graph","digraph")
	 tmp <- charmatch(mode,MODES)
	 if(is.null(tmp)){
	    stop("invalid mode or mode not recognized")
	 }
	 else if(is.na(tmp)){
	    stop("invalid mode or mode not recognized")
	 }
	 else if(tmp==0){
	    stop("ambiguous mode: retry with more characters")
	 }
	 if(degree>=(n-1)) return(Kn(n))
	 if(mode=="graph") 
	    return(rreg(n,degree,scale=scale))
	 MODES <- c("indegree","outdegree")
	 tmp <- charmatch(dmode,MODES)
	 if(is.null(tmp)){
		 stop("invalid dmode or dmode not recognized")
	 }
	 else if(is.na(tmp)){
		 stop("invalid dmode or dmode not recognized")
	 }
	 else if(tmp==0){
		 stop("ambiguous dmode: retry with more characters")
	 }
	 dmode <- MODES[tmp]
	 ed <- c(rep(1,degree),rep(0,n-degree-1))
	 A <- NULL
	 r <- runif(n-1)
	 b <- order(r)
	 z <- c(0,ed[b])
	 if(dmode=="outdegree"){
		 A <- rbind(A,z)
	 }
	 else if(dmode=="indegree"){
		 A <- cbind(A,z)
	 }
	 for(i in 2:(n-1)){
		 r <- runif(n-1)
		 b <- order(r)
		 z <- ed[b]
		 z <- c(z[1:(i-1)],0,z[i:(n-1)])
		 if(dmode=="outdegree"){
			 A <- rbind(A,z)
		 }
		 else if(dmode=="indegree"){
			 A <- cbind(A,z)
		 }
	 }
	 r <- runif(n-1)
	 b <- order(r)
	 z <- c(ed[b],0)
	 if(dmode=="outdegree"){
		 A <- rbind(A,z)
	 }
	 else if(dmode=="indegree"){
		 A <- cbind(A,z)
	 }
	 if(mode=="graph")
		 g <- as.graph(A)
	 else
		 g <- as.digraph(A)
	 g
}

rreg <- function(n,degree=2,scale=T)
{
	 if(((n*degree) %%2) != 0) stop("No such graph possible: nd must be even.")
	 if(degree>=(n-1)) return(Kn(n))
	 NC <- n*(n-1)/2
	 combs <- matrix(0,nrow=NC,ncol=2)
	 k <- 1
	 for(v in 1:(n-1)){
		 for(w in (v+1):n){
			 combs[k,] <- c(v,w)
			 k <- k+1
		 }
	 }
	 p <- rep(degree^2,NC)
	 A <- matrix(0,ncol=n,nrow=n)
	 j <- 1
	 op <- rep(0,NC)
	 while(T){
		 if(!scale) p <- p>0
		 if(all(p==0)){
		    d <- apply(A,1,sum)
			 if(all(d==degree)){
			    break
			 }
			 else{
				 p <- rep(degree^2,NC)
				 A <- matrix(0,ncol=n,nrow=n)
				 j <- j+1
			 }
		 }
		 else if(sum(p==0)==(NC-1)){
		    i <- which(p/sum(p)==1)
			 v <- combs[i,1]
			 w <- combs[i,2]
			 if(A[v,w]==1){
				 p <- rep(degree^2,NC)
				 A <- matrix(0,ncol=n,nrow=n)
				 j <- j+1
			 }
		 }
		 i <- sample(NC,1,prob=p/sum(p))
		 v <- combs[i,1]
		 w <- combs[i,2]
		 A[v,w] <- 1
		 A[w,v] <- 1
		 d <- degree-apply(A,1,sum)
		 for(i in 1:NC){
			 v <- combs[i,1]
			 w <- combs[i,2]
			 p[i] <- d[v]*d[w]
		 }
	 }
	 g <- as.graph(A)
	 g$iterations <- j
	 g
}

rgm <- function(g,m)
{
	n <- graph.order(g)
	mode <- g$mode
	if(mode=="graph" & is.digraph(g)) stop("mode must match graph mode")
	if(mode=="digraph" & is.graph(g)) stop("mode must match graph mode")
	A <- g$A
	if(mode=="graph"){
		A[upper.tri(A)] <- -1
		diag(A) <- -1
		a <- as.vector(A)
		b <- sample((1:length(a))[a==0],m)
		a[b] <- 1
		A <- matrix(a,nrow=n,ncol=n)
		A[A<0] <- 0
		A <- A | t(A)
	}
	else {
	   diag(A) <- -1
		a <- as.vector(A)
		b <- sample((1:length(a))[a==0],m)
		a[b] <- 1
		A <- matrix(a,nrow=n,ncol=n)
		A[A<0] <- 0
	}
	if(mode=="graph")
		g <- as.graph(A)
	else
	   g <- as.digraph(A)
	g
}

rig <- function(n,m,p,withprob=F,threshold=1,tiny=F)
{
   R <- list(0)
	for(i in 1:n){
	    a <- rbinom(m,1,p)
		 if(sum(a)>0) R[[i]] <- (1:m)[a==1]
		 else R[[i]] <- paste("null",i,sep="")
	}
	intersection.graph(R=R,m=m,withprob=withprob,threshold=threshold,tiny=tiny)
}

rstar <- function(n,v,p=1,mode="graph")
{
   g <- empty.graph(n)
	if(missing(v)) v <- sample(1:n,1)
	g$A[v,] <- rbinom(n,1,p)
	if(mode=="graph")
	   g$A[,v] <- g$A[v,]
   else{
	   g$A[,v] <- rbinom(n,1,p)
		g$mode <- "digraph"
	}
	g$A[v,v] <- 0
	g$plot.mode <- "COORDS"
	coord <- cbind(sin(2 * pi * ((0:(n - 2))/(n-1))),
						  cos(2 * pi * ((0:(n - 2))/(n-1))))
	g$coord <- rbind(coord[1:(v-1),],c(0,0),coord[v:(n-1),])
	g
}

rvg <- function(n,r=1,m=1,w=rep(1,m),FUN=NULL)
{
	if(is.null(n)) stop("must provide n")
	if(is.null(FUN))
		vects <- matrix(runif(n*m),ncol=m)
	else
	   vects <- FUN(n,m)
	if(m>1){
		for(i in 1:n){
			v <- sqrt(vects[i,] %*% vects[i,])
			vects[i,] <- vects[i,]/v
		}
	}
	vector.product.graph(vects,r=r,w=w)
}

rtree <- function(n)
{
	if(n<2) stop("n must be at least 2")
	else if(n==2){
	   A <- matrix(1,2,2)
		diag(A) <- 0
		return(as.graph(A))
	}
   sigma <- sample(1:(n-2),replace=T)
	A <- matrix(0,nrow=n,ncol=n)
	S <- 1:n
	for(i in 1:(n-2)){
	   j <- min(setdiff(S,sigma))
		A[j,sigma[1]] <- 1
		sigma <- sigma[2:length(sigma)]
		S <- setdiff(S,j)
	}
	A[S[1],S[2]] <- 1
	g <- as.graph(A)
	g$plot.mode <- "tree"
	g
}

rba <- function(g,m=1,t=10)
{
	pm <- g$plot.mode
	for(i in 1:t){
	   deg <- degrees(g)
		p <- deg/sum(deg)
		n <- graph.order(g)
		g <- add.vertex(g)
		a <- sample(n,m,replace=F,prob=p)
		M <- adjacency.matrix(g)
		M[n+1,a] <- 1
		M[a,n+1] <- 1
		g <- as.graph(M)
	}
	g$plot.mode <- pm
	g
}

rcf <- function(alpha=.5,beta=.5,delta=.5,gamma=.5,p,q,t=10)
{
   g <- empty.graph(1)
	for(i in 1:t){
	   x <- runif(4)
		n <- graph.order(g)
		if(x[1]<alpha && n>1){ # old node
		   if(missing(q)) q <- rep(1,n)
			if(x[3]<delta){
			   v <- sample(n,1)
			}
			else{
				deg <- degrees(g)
				if(sum(deg)==0) deg <- rep(1,n)
				pr <- deg/sum(deg)
			   v <- sample(n,1,prob=pr)
			}
			lq <- min(n,length(q))
			if(sum(q[1:lq])==0) pr <- rep(1/lq,lq)
			else pr <- q[1:lq]/sum(q[1:lq])
			s <- sample(lq,1,prob=pr)
			if(x[4]<gamma){
			   w <- sample(n,s,replace=F)
			}
			else {
				deg <- degrees(g)
				if(sum(deg)==0) deg <- rep(1,n)
				pr <- deg/sum(deg)
			   w <- sample(n,s,replace=F,prob=pr)
			}
			M <- adjacency.matrix(g)
			M[v,w] <- 1
			M[w,v] <- 1
			g <- as.graph(M)
		}
		else{        # new node
			if(missing(p)) p <- rep(1,n)
		   g <- add.vertex(g)
			lp <- min(n,length(p))
			if(sum(p[1:lp])==0) pr <- rep(1/lp,lp)
			else pr <- p[1:lp]/sum(p[1:lp])
			s <- sample(lp,1,prob=pr)
			if(x[2]<beta){
			   v <- sample(n,s,replace=F)
			}
			else{
				deg <- degrees(g)[1:n]
				if(sum(deg)==0) deg <- rep(1,n)
				pr <- deg/sum(deg)
			   v <- sample(n,s,replace=F,prob=pr)
			}
			M <- adjacency.matrix(g)
			M[v,n+1] <- 1
			M[n+1,v] <- 1
			g <- as.graph(M)
		}
	}
	g
}

randomize.graph <- function(g,mode=graph.mode(g),flips=100)
{
    m <- pmatch(mode,c("graph","digraph"))
	if(m==1) g <- as.graph(g)
	else g <- as.digraph(g)
	A <- adjacency.matrix(g)
	n <- graph.order(g)
	v <- sample(n,flips,replace=T)
	w <- sample(n,flips,replace=T)
	for(i in 1:flips){
	   A[v[i],w[i]] <- !A[v[i],w[i]]
	   if(m==1)
		  A[w[i],v[i]] <- A[v[i],w[i]]
	}
	if(m==1) g <- as.graph(A)
	else g <- as.digraph(A)
	g
}
rng <- function(x=NULL,dx=NULL,r=1, p=2,tiny=F, usedeldir=T)
{
   if(is.null(dx)) {
	  if(is.null(x)) stop("One of x or dx must be given.")
	  dx <- pdist(x,p=p)
   }
   n <- nrow(dx)
   A <- matrix(0,nrow=n,ncol=n)
   if(is.vector(x)) x <- matrix(x,ncol=1)
   if(usedeldir && ncol(x)==2 && (length(.find.package("deldir",quiet=T))>0)){
	  require(deldir)
	  del <- deldir(x[,1],x[,2])
	  for(edge in 1:nrow(del$delsgs)){
	     i <- del$delsgs[edge,5]
	     j <- del$delsgs[edge,6]
		 d <- min(apply(cbind(dx[i,-c(i,j)],dx[j,-c(i,j)]),1,max))
		 if(r*dx[i,j] < d){
		    A[i,j] <- 1
		    A[j,i] <- 1
		 }
	  }
   }
   else{
	  diag(dx) <- Inf
	  for(i in 1:n){
		 for(j in setdiff(1:n,i)){
			d <- min(apply(cbind(dx[i,-c(i,j)],dx[j,-c(i,j)]),1,max))
			if(r*dx[i,j] < d){
			   A[i,j] <- 1
			   A[j,i] <- 1
			}
		 }
	  }
   }
	if(tiny){
		out <- as.graph(A)
	}
	else{
		out <- as.graph(A)
		out$plot.mode <- "COORDS"
		out$coord <- x
		out$dx <- dx
		out$r <- r
		out$p <- p
		class(out) <- c("RNG","graph")
	}
   out
}

# Code from Cuppan
rot4d <- function(theta,i,j){
#  Rotation matrix
#       pre-multiplication
#       right hand coordinates  
#       Note: the routine also works for creating a 4-D rotation matrix
#       
	 a <- cos(theta)
	 b <- sin(theta)
	 if(i==2 & j==3)b <- -b
	 mat <- diag(4)
	 mat[i,i] <- a
	 mat[i,j] <- -b
	 mat[j,i] <- b
	 mat[j,j] <- a
	 mat
}

#
# rotates the graph through three dimensions
rotate.graph <- function(g,z=NULL,coord=NULL,
	 theta=0,i=1,j=2,...)
{
	dat <- g$A
	n <- nrow(dat)
	if(is.null(coord))
	   coord <- plot.graph(g,plotit=F,...)
	if(is.null(z)){
	   z <- sqrt(coord[,1]^2+coord[,2]^2)
		z[1:(n/2)] <- -z[1:(n/2)]
	}
   hcoord <- rbind(t(coord),z,rep(1,n))
	plot.graph(g,coord=t(hcoord[1:2,]),...)
	R <- rot4d(theta,i,j)
	p1 <- R %*% hcoord
	plot.graph(g,coord=t(p1[1:2,]/p1[4,]), ...)
}

kneighborhood <- function(g,v=1:graph.order(g),k=1)
{
   if(k<=0) return(v)
   n <- graph.order(g)
   nbhd <- list(0)
	d <- geodesic.dist(g)
   for(i in 1:length(v)){
      nbhd[[i]] <- (1:n)[d[v[i],]<=k]
   }
	if(length(v)==1) nbhd <- unlist(nbhd)
	nbhd
}

hop.disk <- function(g,v=1,hops=1)
{
   if(hops<=0) return(v)
   n <- graph.order(g)
   nbhd <- v
   for(i in 1:hops){
      nbhd <- unique(c(nbhd,neighborhood(g,nbhd,as.indices=T)))
   }
   sort(nbhd)
}

hop.ring <- function(g,v=1,min.hops=1,max.hops=2)
{
   n <- graph.order(g)
   n1 <- hop.disk(g,v,min.hops)
   n2 <- hop.disk(g,n1,max.hops-min.hops)
   setdiff(n2,n1)
}

local.density <- function(g,k=1)
{
	A <- adjacency.matrix(g)
	n <- graph.order(g)
	out <- rep(0,n)
	for(i in 1:n){
	   a <- hop.disk(g,i,hops=k)
		m <- length(a)
		if(m>1)
			out[i] <- sum(A[a,a]>0)/(m*(m-1))
	}
	out
}

local.size <- function(g,k=1)
{
	A <- adjacency.matrix(g)
	n <- graph.order(g)
	out <- rep(0,n)
	for(i in 1:n){
	   a <- hop.disk(g,i,hops=k)
		m <- length(a)
		if(m>1)
			out[i] <- sum(A[a,a]>0)/2
	}
	out
}

local.scan <- function(g,k=1,FUN=graph.size,usehops=T)
{
	if(!usehops)
		a <- kneighborhood(g,k=k)
	n <- graph.order(g)
	out <- rep(0,n)
	for(i in 1:n){
		if(usehops){
		   a <- hop.disk(g,i,hops=k)
			out[i] <- FUN(subgraph(g,a))
		}
		else
			out[i] <- FUN(subgraph(g,a[[i]]))
	}
	out
}


graph.stack <- function(graphs)
{
	if(is.graph(graphs) | is.digraph(graphs)) 
	   return(adjacency.matrix(graphs))

	g <- graphs[[1]]
	if(!inherits(g,"graph")) stop("graphs must be a list of graphs.")
	A <- adjacency.matrix(g)
	n <- length(graphs)
	d <- dim(A)
	gs <- array(0,dim=c(n,nrow(A),ncol(A)))
	for(i in 1:n){
	   g <- graphs[[i]]
		if(!inherits(g,"graph")) stop("graphs must be a list of graphs.")
		A <- adjacency.matrix(g)
		if(any(dim(A) != d)) stop("graphs must all have the same order")
		gs[i,,] <- A
	}
	gs
}

graph.unstack <- function(gs,mode=NULL,thresh=1)
{
   d <- dim(gs)

   if(length(d)==2){
	  if(is.null(mode)){
		 if(is.symmetric(gs)){
			return(as.graph(gs,thresh=thresh))
		 }
		 else{
			return(as.digraph(gs,thresh=thresh))
		 }
	  }
	  else if(mode=="graph"){
	     return(as.graph(gs,thresh=thresh))
	  }
	  else {
	     return(as.digraph(gs,thresh=thresh))
	  }
   }
   else {
	  g <- list(0)
	  for(i in 1:d[1]){
		 if(is.null(mode)){
			if(is.symmetric(gs)){
			   g[[i]] <- as.graph(gs[i,,],thresh=thresh)
			}
			else{
			   g[[i]] <- as.digraph(gs[i,,],thresh=thresh)
			}
		 }
		 else if(mode=="graph"){
			g[[i]] <- as.graph(gs[i,,],thresh=thresh)
		 }
		 else {
			g[[i]] <- as.digraph(gs[i,,],thresh=thresh)
		 }
	  }
   }
   g
}

graph.spectrum <- function(g,mode="laplacian",zapsmall=F)
{
	MODES <- c("laplacian","adjacency","seidel","L")
	tmp <- charmatch(mode,MODES)
	if(is.null(tmp) | is.na(tmp)) stop("invalid mode or mode not recognised")
	else if(tmp==0) stop("ambiguous mode: retry with more characters")
	mode <- MODES[tmp]
	if(mode=="laplacian") A <- laplacian(g)
	else if(mode=="adjacency") A <- adjacency.matrix(g)
	else if(mode=="seidel") A <- seidel(g)
	else if(mode=="L") A <- L(g)
	if(zapsmall) A <- zapsmall(A)
   eigen(A)
}

graph.index <- function(g,mode="adjacency",zapsmall=F)
{
   a <- graph.spectrum(g,mode,zapsmall=zapsmall)
   k <- a$values[1]
   mult <- sum(abs(a$values-k)<10^-6)
   vec <- a$vectors[,1]
   list(index=k,multiplicity=mult,vector=vec)
}

algebraic.connectivity <- function(G,epsilon=10^-10,zapsmall=F)
{
   a <- graph.spectrum(G,mode="laplacian",zapsmall=zapsmall)
	b <- a$values>epsilon
	m <- match(min(a$values[b]),a$values)
	a$value[m-1]
}

algebraic.best.cut <- function(G,zapsmall=F)
{
   n <- graph.order(G)
	best <- algebraic.connectivity(G,zapsmall=zapsmall)
	j <- 0
	for(i in 1:n){
	   a <- algebraic.connectivity(subgraph(G,-i),zapsmall=zapsmall)
		if(a>best){
		   j <- i
			best <- a
		}
	}
	j
}

algebraic.dense <- function(G,zapsmall=F)
{
   x <- NULL
	S<-algebraic.best.cut(G,zapsmall=zapsmall)
	n <- graph.order(G)
	ind <- 1:n
	g <- subgraph(g,-S)
	for(j in 2:(n-1)){
	   y <- algebraic.best.cut(g,zapsmall=zapsmall)
		if(y==0) break
		S[j] <- ind[-S][y]
		g <- subgraph(g,-y)
		cat(S,"\n")
	}
	S
}

fiedler <- function(G,epsilon=10^-10,mode="laplacian",zapsmall=F)
{
   a <- graph.spectrum(G,mode=mode,zapsmall=zapsmall)
	n <- graph.order(G)
	v <- a$vectors[,n:1]
	lambda <- a$values[n:1]
	b <- lambda>epsilon
	m <- match(min(lambda[b]),lambda)
	list(x=v[,m],y=v[,m+1],vectors=v,values=lambda,index=m,lambda2=lambda[m])
}


# Code from Cuppan
#
projr <- function(vd,heye=1.1){
#  Right Eye Projection matrix
#	pre multiplication
#       right hand coordinates, +z toward the viewer
#       d is the viewing distance in inches
#       heye is half the distance between the eyes.
#
	 mat <- vd*diag(4)
	 mat[3,3] <- 0
	 mat[4,3] <- -1
	 #
	 trans1 <- diag(4)
	 trans1[1,4] <- heye
	 #
	 trans2 <- diag(4)
	 trans2[1,4] <- -heye
	 #
	 mat <- trans1%*%mat%*%trans2
	 mat
}

stereo.graph <- function(g,vd=20,z=NULL,coord=NULL,...)
{
	dat <- g$A
	n <- nrow(dat)
	if(is.null(coord))
	   coord <- plot.graph(g,plotit=F,...)
	if(is.null(z)){
	   z <- sqrt(coord[,1]^2+coord[,2]^2)
		z[1:(n/2)] <- -z[1:(n/2)]
	}
   hcoord <- rbind(t(coord),z,rep(1,n))
   proj.rightmat <- projr(vd,1.1)
   proj.leftmat  <- projr(vd,-1.1)

   p1 <- proj.rightmat %*% hcoord
   p2 <- proj.leftmat %*% hcoord
	par(mfrow=c(1,2))
	par(pty="s")
	par(mar=c(0,0,0,0)+.1)
   plot.graph(g,coord=t(p1[1:2,]/p1[4,]),...,axes=F)
	box()
	par(pty="s")
	par(mar=c(0,0,0,0)+.1)
   plot.graph(g,coord=t(p2[1:2,]/p2[4,]),...,axes=F)
	box()
}

# This software developed by David Marchette
# Naval Surface Warfare Center, Dahlgren
# Division, Code B10.  It may be freely used and 
# modified. All warranties, express or implied, are
# disclaimed. 

struct.equiv <- function(graphs)
{
	 if(is.list(graphs)) {
		 if(inherits(graphs,"graph")) graphs <- list(graphs)
		 graphs <- graph.stack(graphs)
	 }
	 else if(!is.array(graphs)){
	    stop("graphs must be a list of graphs, a graph stack, or an adjacency matrix.")
	 }
	 if(length(dim(graphs))==2){
	    n <- nrow(graphs)
		 R <- 1
		 graphs <- array(graphs,c(1,n,n))
	 }
	 else{
	    n <- dim(graphs)[2]
	    R <- dim(graphs)[1]
	 }
	 avail <- 1:n
	 groups <- list(0)
	 g <- 1
	 while(length(avail)>0){
	    i <-  avail[1]
		 groups[[g]] <- i
		 for(j in setdiff(avail,i)){
			 chk <- T
			 for(k in setdiff(1:n,c(i,j))){
				 for(r in 1:R){
					 if(graphs[r,i,k] != graphs[r,j,k]){
						 chk <- F
						 break
					 }
					 if(graphs[r,k,i] != graphs[r,k,j]){
						 chk <- F
						 break
					 }
				 }
				 if(!chk) break
			 }
			 if(chk){
				 groups[[g]] <- c(groups[[g]],j)
				 avail <- setdiff(avail,j)
			 }
		 }
		 avail <- setdiff(avail,i)
		 g <- g+1
	 }
	 groups
}


euclidean.equiv <- function(graph)
{
   A <- adjacency.matrix(graph)
	n <- nrow(A)
	d <- matrix(0,nrow=n,ncol=n)
	for(i in 1:(n-1)){
	   for(j in (i+1):n){
		   k <- setdiff(1:n,c(i,j))
			d[i,j] <- sum((A[k,i]-A[k,j])^2+(A[i,k]-A[j,k])^2)
			d[j,i] <- d[i,j]
		}
	}
	sqrt(d)
}


threshold.graph <- function(w,thresh=1)
{
	 if(is.vector(w)){
	   g <- as.graph(w %o% w,thresh=thresh)
	 }
	 else if(is.matrix(w)){
		if(is.symmetric(w)){
		   g <- as.graph(w,thresh=thresh)
		}
		else{
		   g <- as.digraph(w,thresh=thresh)
		}
	 }
	 else stop("w must be a vector or a matrix")
	 g
}

runif.ball <- function(n,center=c(0,0),radius=1)
{
   thetas <- runif(n,0,2*pi)
   r <- sqrt(runif(n))*radius
   cbind(r*cos(thetas)+center[1],r*sin(thetas)+center[2])
}


symmetrize.graph <- function(g){
	A <- g$A
   g$A <- 1*(A | t(A))
	g
}

is.symmetric <- function(A){
	if(nrow(A) != ncol(A)) return(F)
   !any(A!=t(A))
}

is.graph <- function(g,ignore.symmetry=F){
	if(!inherits(g,"graph")) {
		return(F)
	}
	if(g$mode != "graph") return(F)
	if(ignore.symmetry)
	   return(T)
	else{
		if(is.digraph(g)) return(F)
	   return(is.symmetric(g$A))
	}
}

graph.mode <- function(g)
{
   if(is.digraph(g)) return("digraph")
   else if(is.graph(g)) return("graph")
   else return(NULL)
}

graph.order <- function(g)
{
	A <- adjacency.matrix(g)
	if(nrow(A) != ncol(A)){
	   return(nrow(A)+ncol(A))
	}
	else{
      return(nrow(A))
	}
}

graph.size <- function(g)
{
	A <- adjacency.matrix(g)
	if(nrow(A) != ncol(A)){
	   return(sum(A))
	}
	if(is.digraph(g)){
	   return(sum(A))
	}
	else{
      return(sum(A)/2)
	}
}

graph.density <- function(g)
{
	n <- graph.order(g)
	s <- graph.size(g)
	if(is.digraph(g)){
	   return(s/(n*(n-1)))
	}
	else{
      return(2*s/(n*(n-1)))
	}
}

as.graph <- function(g,thresh=1)
{
	if(is.vector(g)){
		if(length(g)==1){
			g <- matrix(0,nrow=g,ncol=g)
		}
		else{
			g <- matrix(0,nrow=max(g),ncol=max(g))
		}
	}
	if(is.matrix(g)) {
		g <- (g>=thresh)*1
		if(nrow(g) == ncol(g)) 
			diag(g) <- 0
		g <- list(A=g,mode="graph")
		class(g) <- "graph"
		g <- symmetrize.graph(g)
		return(g)
	}
	if(!inherits(g,"graph")) {
		warning("object is not of class graph")
		return(NULL)
	}
	g$A <- g$A>=thresh
	g <- symmetrize.graph(g)
	g$mode <- "graph"
   g
}

as.digraph <- function(g,thresh=1){
	if(is.vector(g)){
		if(length(g)==1){
			g <- matrix(0,nrow=g,ncol=g)
		}
		else{
			g <- matrix(0,nrow=max(g),ncol=max(g))
		}
	}
	if(is.matrix(g)) {
		if(nrow(g)!=ncol(g)) stop("adjacency matrix must be square for digraphs")
		g <- (g>=thresh)*1
		diag(g) <- 0
		g <- list(A=g,mode="digraph")
		class(g) <- "graph"
		return(g)
	}
	if(!inherits(g,"graph")) {
		warning("object is not of class graph")
		return(NULL)
	}
	g$A <- g$A>=thresh
	g$mode <- "digraph"
   g
}

is.digraph <- function(g,strict=F){
	if(!inherits(g,"graph")) {
		return(F)
	}
	if(g$mode != "digraph") return(F)
	if(strict)
		return(!is.symmetric(g$A))
   else
	   return(T)
}

is.complete <- function(g){
	if(is.matrix(g)) A <- g
	else A <- g$A
   n <- nrow(A)
   sum(A) == n*(n-1)
}

degrees <- function(g,cmode="outdegree")
{
   MODES <- c("indegree","outdegree","total")
	tmp <- charmatch(cmode,MODES)
	if(is.null(tmp) | is.na(tmp)){
	   stop("invalid cmode or cmode not recognized")
	}
	if(tmp==0){
	   stop("ambiguous cmode: retry with more characters")
	}
	mode <- MODES[tmp]
	if(mode=="indegree"){
	   d <- apply(g$A,2,sum)
	}
	else if(mode=="outdegree"){
	   d <- apply(g$A,1,sum)
	}
	else{
	   d <- apply(g$A,1,sum)+apply(g$A,2,sum)
	}
	d
}

degree.sequence <- function(G,cmode="outdegree")
{
   sort(degrees(G,cmode),decreasing=T)
}

# only for labeled graphs! NOT an isomorphism check!
graph.equal <- function(g1,g2)
{
   if(graph.order(g1) != graph.order(g2)) return(FALSE)
   if(graph.size(g1) != graph.size(g2)) return(FALSE)
   sum(g1$A!=g2$A)==0
}

# if it returns TRUE then the graphs are isomorphic
# if it returns FALSE then the graphs are not isomorphic
# if it a NULL then it couldn't tell one way or another
graph.isomorphic.random <- function(g1,g2,tries=100, check.clique=T)
{
	.iso.reason <- ""
	.iso.reorder <<- NULL
	if(graph.order(g1) != graph.order(g2)) {
	   .iso.reason <<- "order"
	   return(FALSE)
	}
	if(graph.size(g1) != graph.size(g2)) {
	   .iso.reason <<- "size"
	   return(FALSE)
	}
	if(number.components(g1) != number.components(g2)) {
	   .iso.reason <<- "number components"
	   return(FALSE)
	}
	x <- component.sequence(g1)
	y <- component.sequence(g2)
	if(any(x!=y)) {
	   .iso.reason <<- "component sequence"
	   return(FALSE)
	}
	d1 <- sort(degrees(g1))
	d2 <- sort(degrees(g2))
	if(any(d1!=d2)) {
	   .iso.reason <<- "degree sequence"
	   return(FALSE)
	}
	if(graph.equal(g1,g2)) {
	   .iso.reorder <<- 1:graph.order(g1)
	   .iso.reason <<- "degree reorder"
	   return(TRUE)
	}
	gorder <- order(degrees(g1))
	g1 <- reorder.vertices(g1)
	g2 <- reorder.vertices(g2)
	if(graph.equal(g1,g2)) {
	   .iso.reorder <<- gorder[order(degrees(g2))]
	   .iso.reason <<- "random"
	   return(TRUE)
	}
	n <- graph.order(g1)
	x <- graph.spectrum(g1)$values
	y <- graph.spectrum(g2)$values
	if(any(abs(x-y)>10^-6)){
	   .iso.reason <<- "spectrum"
	   return(FALSE)
	}
	i <- 1
    while(i <= tries){
	   i <- i+1
	   dd <- d1+runif(n,0,.1)
	   od <- order(dd)
	   g <- reorder.vertices(g1,ord=od)
	   if(graph.equal(g,g2)) {
		  .iso.reorder <<- gorder[od]
		  .iso.reason <<- "random"
		  return(TRUE)
	   }
	}
	if(check.clique){
	   if(clique.number(g1) != clique.number(g2)) {
		  .iso.reason <<- "clique number"
		  return(FALSE)
	   }
	}
	return(NULL)
}

graph.isomorphic <- function(g1,g2)
{
   .iso.reason <<- ""
   .iso.reorder <<- NULL
	q <- graph.isomorphic.random(g1,g2,100,check.clique=T)
	k <- 0
	if(is.null(q)){
	   n <- graph.order(g1)
	   gorder <- order(degrees(g1))
	   g1 <- reorder.vertices(g1)
	   g2 <- reorder.vertices(g2)
	   d <- degrees(g1)
	   num <- (n-1):0
	   while(any(num!=0)){
		  k <- k+1
		  num[1] <- num[1]+1
		  for(i in 1:(n-1)){
		     if(num[i]!=n) break
			 num[i] <- 0
			 num[i+1] <- num[i+1]+1
		  }
		  if(num[n]==n) num[n] <- 0
		  if(length(unique(num))==n){
			 s <- rev(num+1)
			 if(all(d==d[s])){
				q <- graph.equal(g1,graph.reorder(g2,s))
				if(q){
				   .iso.reorder <<- gorder[s]
				   .iso.reason <<- "systematic"
				   return(q)
				}
			 }
			 else {
				dd <- d+runif(n,0,.1)
				s <- order(dd)
				q <- graph.equal(g1,graph.reorder(g2,s))
				if(q){
				   .iso.reorder <<- gorder[s]
				   .iso.reason <<- "random"
				   return(q)
				}
			 }
		  }
		  else {
			 dd <- d+runif(n,0,.1)
			 s <- order(dd)
			 q <- graph.equal(g1,graph.reorder(g2,s))
			 if(q){
				.iso.reorder <<- gorder[s]
				.iso.reason <<- "random"
				return(q)
			 }
		  }
		  if((k %% 10000)==0) cat(k," ")
	   }
   }
   else return(q)
   .iso.reason <<- "exhaustive"
   return(FALSE)
}

"==.graph" <- function(x,y)  graph.isomorphic(x,y)
"!=.graph" <- function(x,y) !graph.isomorphic(x,y)

adjacency.matrix <- function(g,augmented=F){
   A <- g$A
	if(augmented) diag(A) <- 1
	A
}

# a vector indicating which vertices
# have indegree 0
indegree0 <- function(g){
   degrees(g,cmode="indegree")==0
}

reorder.vertices <- function(g,ord=order(degrees(g)))
{
   g$A <- adjacency.matrix(g)[ord,ord]
   if(!is.null(g$coord)) {
	  if(nrow(g$coord)==length(ord))
		 g$coord <- g$coord[ord,]
   }
   if(!is.null(g$x)) {
	  if(is.vector(g$x)) {
		 if(length(g$x)==length(ord))
			g$x <- g$x[ord]
	  }
	  else {
		 if(nrow(x)==length(ord))
			g$x <- g$x[ord,]
	  }
   }
   if(!is.null(g$R)) {
	  if(length(g$R)==length(ord))
		 g$R <- g$R[ord]
   }
   g
}

graph.reorder <- function(g,ord=sample(graph.order(g)))
{
   reorder.vertices(g,ord)
}

# a vector indicating the neighbors
# of a vertex (or set of vertices) v 
neighborhood <- function(g,v=1,open=F,as.indices=F)
{
	A <- adjacency.matrix(g)
   n <- nrow(A)
   nbhd <- rep(0,n)
   nbhd[v] <- !open
   for(i in v)
      nbhd <- nbhd | A[i,]
   if(as.indices) nbhd <- (1:n)[nbhd]
   nbhd
}

is.isolated <- function(g,vertex=1)
{
	all(g$A[vertex,]==0)&all(g$A[,vertex]==0)
}

isolated.vertices <- function(g)
{
	degrees(g,cmode="total")==0
}

remove.isolates <- function(g)
{
	z <- !isolated.vertices(g)
	g$A <- g$A[z,z]
	if(is.graph(g)) g <- as.graph(g$A)
	else if(is.digraph(g)) g <- as.digraph(g$A)
	else stop("g must be a graph or digraph")
	g
}

remove.vertices <- function(g,v=1)
{
	g$A <- g$A[-v,-v]
	if(is.graph(g)) g <- as.graph(g$A)
	else if(is.digraph(g)) g <- as.digraph(g$A)
	else stop("g must be a graph or digraph")
	g
}

bidirected <- function(g)
{
   g$A <- 1*(g$A & t(g$A))
	as.graph(g)
}

remove.bidirected <- function(g)
{
   g$A <- g$A-1*(g$A & t(g$A))
	as.digraph(g)
}

summary.graph <- function(object,verbose=F,...)
{
	g <- object
	if(is.digraph(g)) type <- "digraph"
	else if(is.graph(g)) type <- "graph"
	else stop("g must be a graph")
	n <- graph.order(g)
	d <- degrees(g,cmode="outdegree")
	dg <- dominate(g)
	comp <- graph.components(g)
	num <- matrix(0,ncol=max(comp),nrow=1)
	for(i in 1:max(comp)) {
		num[i] <- sum(comp==i)
	}
	colnames(num) <- paste("Comp",1:max(comp),sep="")
   out <- list(
		  type=type,
	     order=n,
	     size=graph.size(g),
		  density=graph.density(g),
		  gamma.hat=length(dg),
		  dominating.set=dg,
		  num.components=num)
	if(type=="graph") {
		out$degrees=summary(d)
		out$num.isolates=sum(d==0)
	}
	else{
	   out$out.degrees=summary(d)
		out$in.degrees=summary(degrees(g,cmode="indegree"))
		out$num.isolates=sum(num==1)
	}
	if(verbose | n<100){
		if(type=="graph")
			out$clique.number=clique.number(g)
		geo <- geodesic.dist(g)
		out$geodesics <- summary(as.vector(geo))
	}
   out
}

str.graph <- function(object,show.adj=T,...)
{
	if(show.adj){
		A <- adjacency.matrix(object)
		b <- apply(A,1,paste,collapse=" ")
		for(i in 1:length(b))
		   cat("\t",b[i],"\n")
	}
	else str.default(object)
}

delete.vertex <- function(g, v=1)
{
	remove.vertices(g,v)
}

incidence.matrix <- function(g)
{
   e <- graph.edges(g)
   n <- graph.order(g)
   A <- matrix(0,nrow=n,ncol=nrow(e))
   for(i in 1:nrow(e)){
      A[e[i,],i] <- 1
   }
   A
}

graph.edges <- function(g)
{
  A <- adjacency.matrix(g)
  edgelist <- NULL
  n <- nrow(A)
  if(is.digraph(g)){
     for(i in 1:n){
	     a <- (1:n)[A[i,]==1]
		  if(length(a)>0) edgelist <- rbind(edgelist,cbind(rep(i,length(a)),a))
	  }
  }
  else{
     for(i in 2:n){
	     a <- (1:i)[A[i,(1:i)]==1]
		  if(!is.na(a) && length(a)>0) 
		      edgelist <- rbind(edgelist,cbind(rep(i,length(a)),a))
	  }
  }
  if(!is.matrix(edgelist)) edgelist <- matrix(edgelist,ncol=2)
  edgelist
}

delete.edge <- function(g, v=NULL, w)
{
	if(graph.size(g)==0) return(g)
	if(is.null(v)){
	   a <- graph.edges(g)
		n <- nrow(a)
		x <- a[sample(1:n,1),]
		v <- x[1]
		w <- x[2]
	}
	if(is.matrix(v)) {
		w <- v[,2]
		v <- v[,1]
	}
	if(is.graph(g)){
		g$A[cbind(w,v)] <- 0
	}
   g$A[cbind(v,w)] <- 0
	g
}

flip.edge <- function(g,v=NULL, w)
{
	if(is.null(v)){
		n <- graph.order(g)
		x <- sample(1:n,2,replace=F)
		v <- x[1]
		w <- x[2]
	}
	if(is.graph(g)){
	   g$A[cbind(v,w)] <- !g$A[cbind(v,w)]
	   g$A[cbind(w,v)] <- !g$A[cbind(w,v)]
	}
	else{
	   g$A[cbind(v,w)] <- !g$A[cbind(v,w)]
	}
	g
}

add.vertex <- function(g)
{
   M <- adjacency.matrix(g)
	n <- nrow(M)
	M <- cbind(rbind(M,rep(0,n)),rep(0,n+1))
	if(is.digraph(g)) g <- as.digraph(M)
	else g <- as.graph(M)
	g
}

add.edge <- function(g, v=NULL, w)
{
	if(is.null(v)){
	   a <- degrees(g)
		n <- graph.order(g)
		b <- (1:n)[a<n]
		v <- sample(b,1)
		A <- adjacency.matrix(g)
		b <- (1:n)[A[v,]==0]
		w <- sample(b,1)
	}
	if(is.matrix(v)) {
		w <- v[,2]
		v <- v[,1]
	}
	if(is.graph(g)){
		g$A[cbind(w,v)] <- 1
	}
   g$A[cbind(v,w)] <- 1
	g
}

graph.diff <- function(g,h)
{
   if(graph.order(g) != graph.order(h)) 
	   stop("graphs must have the same number of vertices")
	if(graph.size(h)==0) return(g)
   delete.edge(g,graph.edges(h))
}

graph.add <- function(g,h)
{
   if(graph.order(g) != graph.order(h)) 
	   stop("graphs must have the same number of vertices")
   as.graph(g$A | h$A)
}

"+.graph" <- function(x,y) graph.add(x,y)
"-.graph" <- function(x,y) graph.diff(x,y)

"^.graph" <- function(x,n)
{
   if(is.infinite(n)){
      y <- component.sequence(x)
	  if(length(y)==1) return(Kn(y))
	  g <- Kn(y[1])
	  for(i in 2:length(y)){
	     g <- graph.adjoin(g,Kn(y[i]))
	  }
	  return(g)
   }
   y <- as.integer(n)
   if(y==0){
      x$A <- matrix(0,nrow=graph.order(x),ncol=graph.order(x))
	  return(x)
   }
   if(y<=0) stop("power must be a positive integer")
   if(y==1) return(g)
   A <- x$A
   for(i in 2:y)
	  x$A <- 1*((x$A %*% A)>0)
   x
}

graph.adjoin <- function(g,h,v=NULL,mode="graph")
{
	n <- graph.order(g)
	m <- graph.order(h)
   A <- matrix(0,nrow=n+m,ncol=n+m)
	A[1:n,1:n] <-adjacency.matrix(g)
	A[(n+1):(n+m),(n+1):(n+m)] <- adjacency.matrix(h)
	if(!is.null(v)){
		if(is.vector(v)) v <- matrix(v,nrow=1)
		for(i in 1:nrow(v)){
			A[v[i,1],v[i,2]+n] <- 1
		}
	}
	if(mode=="graph") 
		g <- as.graph(A)
	else 
	   g <- as.digraph(A)
   g
}

graph.complement <- function(g)
{
   g$A <- matrix(as.integer(!g$A),nrow=nrow(g$A))
	diag(g$A) <- 0
	g
}

digraph.reverse <- function(g)
{
   g$A <- t(g$A)
	g
}

line.graph <- function(g)
{
   v <- graph.edges(g)
	intersection.graph(lapply(1:nrow(v),function(i)v[i,]))
}

#prisner
gallai.graph <- function(g)
{
   v <- graph.edges(g)
	h <- line.graph(g)
	A <- adjacency.matrix(h)
	B <- adjacency.matrix(g)
	n <- nrow(v)
	for(i in 1:(n-1)){
	   for(j in (i+1):n){
			k <- unique(as.vector(v[c(i,j),]))
			if((B[k[1],k[2]]==1) && (B[k[1],k[3]]==1) && (B[k[2],k[3]]==1)){
				A[i,j] <- 0
				A[j,i] <- 0
			}
		}
	}
	as.graph(A)
}

product.graph <- function(g,h)
{
   n <- graph.order(g)
   m <- graph.order(h)
	G <- adjacency.matrix(g)
	H <- adjacency.matrix(h)
	A <- matrix(0,nrow=n*m,ncol=n*m)
	for(i in 1:n){
		for(j in 1:m){
			for(k in 1:n){
				for(l in 1:m){
					if((i==k && H[j,l]) || (j==l && G[i,k])){
					   A[(i-1)*m+j,(k-1)*m+l] <- 1
					}
				}
			}
		}
	}
	diag(A) <- 0
	if(is.graph(g) && is.graph(h)) g <- as.graph(A)
	else g <- as.digraph(A)
	g
}

is.bipartite <- function(g)
{
	 n <- graph.order(g)
	 a <- 1
	 G <- symmetrize.graph(g)
	 b <- neighborhood(G,a,as.indices=T,open=T)
	 while(T){
		 a <- unique(neighborhood(G,b,as.indices=T,open=T))
		 b <- unique(neighborhood(G,a,as.indices=T,open=T))
		 if(length(c(a,b))>=n) break
	 }
	 length(intersect(a,b))==0
}

graph.coords <- function(g)
{
	coords <- g$coords
   if(is.null(coords)) coords <- plot(g,plotit=F)
	coords
}

is.eulerian <- function(g)
{
   is.graph(g) && is.connected(g) && all((degrees(g) %% 2)==0)
}

graph.join <- function(g,h)
{
	grph <- is.graph(g) & is.graph(h)
	n <- graph.order(g)
	m <- graph.order(h)
   k <- graph.adjoin(g,h,NULL)
	A <- adjacency.matrix(k)
	B <- matrix(1,nrow=n,ncol=m)
	A[1:n,(n+1):(n+m)] <- 1
	A[(n+1):(n+m),1:n] <- 1
	if(grph)
		return(as.graph(A))
	else
		return(as.digraph(A))
}

subgraph <- function(g,v)
{
	A <- adjacency.matrix(g)
	h <- g
	h$A <- matrix(A[v,v],nrow=length(v),ncol=length(v))
	if(!is.null(g$coord)) h$coord <- g$coord[v,]
	h
}

outdegree.distribution <- function(g,plotit=T,
                                   xlab="Outdegree",ylab="Number of Vertices",
											  ...)
{
   deg <- degrees(g,cmode="outdegree")
   out <- table(deg)
	x <- sort(unique(deg))
	y <- as.vector(table(deg))
	if(plotit){
	   plot(x,y,log="xy",xlab=xlab,ylab=ylab,...)
	}
	invisible(list(x=x,y=y,degrees=deg))
}

cyclomatic.number <- function(g)
{
   graph.size(g)-graph.order(g)+1
}
.First.lib<- function(lib,pkg) {
	library.dynam("graph",pkg,lib)
}
require(stats)
cat("graph utilities\nVersion: 0.3-16.7\n\n")
if(!exists(".graph.loaded.once")){
   cat("For full capabilities, get R packages: ")
   cat("deldir, gregmisc, sna\n")
   cat("and programs qhull (www.qhull.org) and\n")
   cat("graphviz (www.research.att.com/sw/tools/graphviz)\n")
   .graph.loaded.once <<- T
}

mpoly <- function(classes)
{
   u <- unique(classes)
   coord <- matrix(0,nrow=length(classes),ncol=2)
   nc <- length(u)
   theta <- seq(0,2*pi,length=nc+1)
   X <- sin(theta)
   Y <- cos(theta)
   if(nc==2) {
     coord=matrix(data=0,nrow=length(classes),ncol=2)
     coord[,1]=classes
     m <- sum(classes==u[1])
     coord[classes==u[1],2]=(1:m)
     m <- sum(classes==u[2])
     coord[classes==u[2],2]=(1:m)
   } else{
      for(i in 1:nc){
         m <- sum(classes==u[i])
         xc <- as.matrix(seq(X[i],X[i+1],length=m+2)[1:m+1])
         yc <- as.matrix(seq(Y[i],Y[i+1],length=m+2)[1:m+1])
         zc <- cbind(xc,yc)
         coord[classes==u[i],] <- zc
      }
   }
   coord
}
mflower <- function(classes, R=length(unique(classes)),r=1,
				center=NULL)
{
   u <- unique(classes)
   coord <- matrix(0,nrow=length(classes),ncol=2)
   nc <- length(u)
	if(length(r)==1) r <- rep(r,nc)
	if(is.null(center)) {
		center <- -10^10
		theta <- seq(0,2*pi,length=nc+1)
	}
	else{
		if(any(classes==center))
			theta <- seq(0,2*pi,length=nc)
		else
			theta <- seq(0,2*pi,length=nc+1)
	}
	j <- 1
   for(i in 1:nc){
      m <- sum(classes==u[i])
      phi <- seq(theta[i]-pi/2,theta[i]+pi/2,length=m+1)[1:m]
      cr <- R+r[i]
      x1 <- c(cr*cos(theta[j]),cr*sin(theta[j]))
		if(u[i]==center){
			coord[classes==u[i],] <- cbind(r[i]*cos(phi),r[i]*sin(phi))
		}
		else {
			coord[classes==u[i],] <- cbind(r[i]*cos(phi)+x1[1],r[i]*sin(phi)+x1[2])
			j <- j+1
		}
   }
	coord
}

nicelabel <- function (x, gmode = NULL, 
    label,label.coord=NULL,
    coord = NULL, jitter = FALSE, usearrows = TRUE, 
    mode = "circle", displayisolates = TRUE, pad = 0, vertex.pch = 19, 
    label.cex = 1, vertex.cex = 1, label.col = 1, edge.col = 1, 
    vertex.col = 1, arrowhead.length = 0.2, edge.type = 1, edge.lwd = 0, 
	 xlab=NULL,ylab=NULL,
	 newplot=T, arrowtype="arrow",head.lty=edge.type,head.lwd=edge.lwd,
	 head.col=edge.col+1,xlim,ylim,plotit=T,m=2,rot=F,
	 clique.cex=2,clique.label.cex=label.cex[1],clique.label.col=label.col[1],
	 keep.circle=F,clique.num=NULL,shadow=F,
	 cspring=1,crep=1,l=1,iterations=10,delta=0.01,start="random",
	 parts=NULL,force.mode=F,
	 bycomponent=F,component=0, 
	 R=3, r=1,center=NULL,
	 foffset=0,
	 ...)
{
	 g <- x
	 if(bycomponent){
	    g <- graph.component(g,component)
	 }
	 if(force.mode) g$plot.mode <- NULL
	 if(!force.mode & !is.null(g$plot.mode) & is.null(coord)){
		 mode <- g$plot.mode
		 if(!is.null(g$xlim)) xlim <- g$xlim
		 if(!is.null(g$ylim)) ylim <- g$ylim
		 if(mode=="bipartite") parts <- g$parts
		 if(mode=="pinwheel") parts <- g$parts
		 if(mode=="mpoly") parts <- g$parts
		 if(mode=="mflower") parts <- g$parts
		 else if(mode=="mcircles"){
			 m <- g$m
			 if(!is.null(g$rot)) rot <- g$rot
		 }
		 else if(mode=="COORDS"){
			 coord <- g$coord
			 if(is.null(coord))
			    mode <- "circle"
			 else
				 mode <- NULL
			 if(is.vector(coord)){
			    coord <- cbind(coord,runif(length(coord)))
			 }
		 }
	 }
	 if(!is.graph(g,ignore.symmetry=T) && !is.digraph(g)){
		 stop("object is not a graph")
	 }
	 if(missing(label)) label <- 1:graph.order(g)
	 ARROWS = c("arrow","midarrow","endarrow")
	 arrowtype <- pmatch(arrowtype,ARROWS)
	 if(arrowtype<1 || arrowtype>3){
	    warning("unknown arrowtype. using type \"arrow\"")
		 arrowtype<-1
	 }
	 d <- g$A
	 if(is.null(gmode)){
	    if(is.graph(g)) gmode <- "graph"
		 else gmode <- "digraph"
	 }
    if (gmode == "graph") {
        usearrows <- FALSE
    }
    n <- nrow(d)
    d[is.na(d)] <- 0
	 if(!is.null(mode)){
		 MODES <- c("princoord","eigen","seidel","laplacian","fiedler","mds","random","circle",
						"circrand","rmds","geodist","mst","dot","adj","seham","mcircles",
						"spring","tree","clique","bipartite","line","pinwheel","mpoly","mflower")
		 tmp <- charmatch(mode,MODES)
		 if(is.null(tmp)){
			 stop("invalid mode or mode not recognized")
		 }
		 else if(is.na(tmp)){
			 stop("invalid mode or mode not recognized")
		 }
		 else if(tmp==0){
			 stop("ambiguous mode: retry with more characters")
		 }
		 else if(tmp<1 || tmp>length(MODES)){
			 stop("invalid mode or mode not recognized")
		 }
		 mode <- MODES[tmp]
		 if(mode %in% 
			 c("princoord","eigen","seidel","laplacian","fiedler","adj","mds","sehan")){
				 if(is.null(xlab))
					 xlab <- expression(lambda[1])
				 if(is.null(ylab))
					 ylab <- expression(lambda[2])
			 }
		 else{
			 if(is.null(xlab))
				 xlab <- ""
			 if(is.null(ylab))
				 ylab <- ""
		 }
	 }
    if (!is.null(coord)) {
        x <- coord[, 1]
        y <- coord[, 2]
    }
	 else if (mode == "pinwheel"){
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
	    coord <- pinwheel(classes=parts,R=R,r=r,center=center)
            
		 x <- coord[,1]
		 y <- coord[,2]
	 }
	 else if (mode == "mpoly"){
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
	    coord <- mpoly(classes=parts)
		 x <- coord[,1]
		 y <- coord[,2]
	 }
	 else if (mode == "mflower"){
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
	    coord <- mflower(classes=parts)
		 x <- coord[,1]
		 y <- coord[,2]
	 }
    else if (mode == "princoord") {
        cd <- cor(rbind(d, t(d)), use = "pairwise.complete.obs")
        cd <- replace(cd, is.na(cd), 0)
        e <- eigen(cd, symmetric = TRUE)
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "eigen") {
        e <- eigen(d, symmetric = (gmode != "digraph"))
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "laplacian") {
        e <- eigen(laplacian(g), symmetric = (gmode != "digraph"))
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "fiedler") {
        e <- fiedler(g)
        x <- Re(e$vectors[, e$index+foffset])
        y <- Re(e$vectors[, e$index+1+foffset])
    }
    else if (mode == "seidel") {
        e <- eigen(seidel(g), symmetric = (gmode != "digraph"))
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "mds") {
        Dmat <- matrix(nrow = n, ncol = n)
        diag(Dmat) <- 0
        for (i in 1:n) for (j in 1:n) if (i > j) 
            Dmat[i, j] <- sqrt(sum(abs(d[i, ] - d[j, ])) + sum(abs(d[, 
                i] - d[, j])))
        Dmat[upper.tri(Dmat)] <- t(Dmat)[upper.tri(Dmat)]
        Imat <- matrix(nrow = n, ncol = n)
        Imat[, ] <- 0
        diag(Imat) <- 1
        Hmat <- Imat - (1/n) * rep(1, n) %*% t(rep(1, n))
        Amat <- -0.5 * Dmat^2
        Bmat <- Hmat %*% Amat %*% Hmat
        e <- eigen(Bmat)
        x <- Re(e$vectors[, 1])
        y <- Re(e$vectors[, 2])
    }
    else if (mode == "random") {
        x <- runif(n, -1, 1)
        y <- runif(n, -1, 1)
    }
    else if (mode == "line") {
        x <- 1:n
        y <- rep(0,n)
    }
    else if (mode == "circle") {
        x <- sin(2 * pi * ((0:(n - 1))/n))
        y <- cos(2 * pi * ((0:(n - 1))/n))
    }
    else if (mode == "circrand") {
        tempd <- rnorm(n, 1, 0.25)
        tempa <- runif(n, 0, 2 * pi)
        x <- tempd * sin(tempa)
        y <- tempd * cos(tempa)
    }
    else if (mode == "rmds") {
        require(mva)
        tempmds <- cmdscale(dist(d))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "geodist") {
        require(mva)
        tempmds <- cmdscale(as.dist(geodesic.dist(g)))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "mst") {
		  q <- geodesic.dist(g)
        temp <- mst(distmat=q)
		  z <- dotlayout(temp)
        x <- z[, 1]
        y <- z[, 2]
    }
    else if (mode == "dot") {
		  z <- dotlayout(g)
        x <- z[, 1]
        y <- z[, 2]
    }
    else if (mode == "adj") {
        require(mva)
        tempmds <- cmdscale(as.dist(-d + max(d)))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "seham") {
        require(mva)
        require(sna)
        temp <- sedist(d)
        tempmds <- cmdscale(as.dist(temp))
        x <- tempmds[, 1]
        y <- tempmds[, 2]
    }
    else if (mode == "mcircles") {
		  if(m<2 || m>(n/2)) stop("invalid choice of m in mcircles")
		  r <- m:1
		  x <- NULL
		  y <- NULL
		  N <- n/m
		  for(i in 1:m){
			  x <- c(x,r[i]*sin(2 * pi * ((0:(N - 1))/N)-rot*(i %% 2)*pi/N))
			  y <- c(y,r[i]*cos(2 * pi * ((0:(N - 1))/N)-rot*(i %% 2)*pi/N))
		  }
		  N <- n-length(x)
		  if(N>0){
			  x <- c(x,r[m]*sin(2 * pi * ((0:(N - 1))/N)-rot*((m+1) %% 2)*pi/N))
			  y <- c(y,r[i]*cos(2 * pi * ((0:(N - 1))/N)-rot*((m+1) %% 2)*pi/N))
		  }
    }
	 else if (mode=="spring"){
		if(is.matrix(start)){
			if(nrow(start) != n || ncol(start) != 2) stop("start must be nx2")
			points <- start
		}
		else if(is.character(start)){
			points <- Recall(g, gmode = gmode, mode="random",plotit=F)
	 }
	 else stop("start must be either a mode name or a matrix")
		frep <- function(x,y,crep){
			dd <- pdist(x,y,d=2)
			crep/dd^3 * (x-y)
		}
		fspring <- function(x,y,cspring,l){
			dd <- pdist(x,y,d=2)
			cspring*log(dd/l)*(x-y)
		}
		force <- function(A,j,points,crep,cspring,l){
			f <- 0
			for(i in setdiff(1:nrow(A),j)){
				if(A[i,j]==1){
					f <- f+fspring(points[i,],points[j,],cspring,l)
				}
				else{
					f <- f+frep(points[i,],points[j,],crep)
				}
			}
			f
		}
		direction <- function(x,y){
			#acos(x %*% y/(sqrt(sum(x %*% x)*sum(y %*% y))))
			atan2(x ,y)
		}
		B <- d
		diag(B) <- 0
		oforces <- matrix(0,nrow=n,ncol=2)
		deltas <- rep(0,n)
		for(i in 1:iterations){
			forces <- matrix(0,nrow=n,ncol=2)
			for(j in 1:n){
				forces[j,] <- force(B,j,points,crep,cspring,l)
				if(i>1)
					dd <- direction(forces[j,],oforces[j,])
				else 
					dd <- 1
				if((dd>0) && (dd<.01)) {
					deltas[j] <- delta*2
				}
				else if((dd<0) && (dd> -0.01)) {
					deltas[j] <- delta/2
				}
				else {
					deltas[j] <- delta
				}
				oforces[j,] <- forces[j,]
			}
			for(j in 1:n){
				points[j,] <- points[j,]+delta*forces[j,]
			}
			#plot.graph(d,coord=points)
		}
		x <- points[,1]
		y <- points[,2]
	 }
	 else if(mode=="tree"){
	    v <- 1
		 x <- rep(0,n)
		 y <- rep(0,n)
		 S <- 1:n
		 i <- 1
		 while(length(S)>0){
			 q <- length(v)
			 k <- -q
			 for(j in 1:length(v)){
				 a <- neighborhood(g,v[j],open=T,as.indices=T)
				 b <- setdiff(a,v)
				 if(length(b)>0){
					 S <- setdiff(S,a)
					 z <- seq(k,k+1,length=length(b))
					 v <- c(v,a)
					 x[b] <- z
					 y[b] <- i
				 }
				 k <- k + length(b)+1
			 }
			 i <- i+1
		 }
	 }
	 else if(mode=="clique"){
		 if(is.complete(g)){
		    stop("graph is complete")
		 }
		 a <- clique.number(g)
		 if(!is.null(clique.num)){
			 if(a<clique.num) 
			    stop("clique.num must be less than the clique number of the graph")
		    a <- clique.num
		 }
		 if(a < n && a>2){
			 g <- collapse(g,vertices=clique(g,a))
			 if(!is.null(label)){
			    if(shadow){
			       slabel <- label[g$collapsed]
				 }
			    label <- c(paste(label[setdiff(1:n,g$collapsed)]),
				            eval(substitute(
								     expression(K[k]),list(k=length(g$collapsed)))))
			 }
			 d <- g$A
			 N <- nrow(g$A)
			 ov <- vertex.cex
			 vertex.cex <- rep(ov,N)
			 vertex.cex[N] <- clique.cex+ov
			 if(length(vertex.pch)>1 & length(vertex.pch)<n){
			    vertex.pch <- rep(vertex.pch,n/length(vertex.pch) +1)[1:n]
			 }
			 if(length(vertex.pch)==n){
			    vertex.pch <- vertex.pch[setdiff(1:n,g$collapsed)]
				 vertex.pch <- c(vertex.pch,13)
			 }
			 else{
				 vertex.pch <- rep(vertex.pch,N)
				 vertex.pch[N] <- 13
			 }
			 if(length(label.col)>1 & length(label.col)<n){
			    label.col <- rep(label.col,n/length(label.col) +1)[1:n]
			 }
			 if(length(label.col)==n){
			    label.col <- label.col[setdiff(1:n,g$collapsed)]
				 label.col <- c(label.col,clique.label.col)
			 }
			 else{
				 label.col <- rep(label.col,N)
				 label.col[N] <- clique.label.col
			 }
			 if(length(label.cex)>1 & length(label.cex)<n){
			    label.cex <- rep(label.cex,n/length(label.cex) +1)[1:n]
			 }
			 if(length(label.cex)==n){
			    label.cex <- label.cex[setdiff(1:n,g$collapsed)]
				 label.cex <- c(label.cex,clique.label.cex)
			 }
			 else{
				 label.cex <- rep(label.cex,N)
				 label.cex[N] <- clique.label.cex
			 }
			 if(keep.circle){
				  a <- sin(2 * pi * ((0:(n - 1))/n))
				  b <- cos(2 * pi * ((0:(n - 1))/n))
				  x <- c(a[setdiff(1:n,g$collapsed)],0)
				  y <- c(b[setdiff(1:n,g$collapsed)],0)
				  if(shadow){
				     z <- a[g$collapsed]
				     w <- b[g$collapsed]
				  }
			 }
			 else{
				 m <- N-1
				 x <- sin(2 * pi * ((0:(m - 1))/m))
				 y <- cos(2 * pi * ((0:(m - 1))/m))
				 x <- c(x,0)
				 y <- c(y,0)
			 }
			 if(!is.null(label.coord)){
				 xl <- c(label.coord[setdiff(1:n,g$collapsed),1],0)
				 yl <- c(label.coord[setdiff(1:n,g$collapsed),2],0)
				 if(shadow){
					 z <- c(label.coord[g$collapsed,1],0)
					 w <- c(label.coord[g$collapsed,2],0)
				 }
			 }
			 n <- N
		 }
		 else {
		    stop("Either g is complete, or contains no cliques")
		 }
	 }
	 else if(mode=="bipartite"){
		 if(is.null(parts)){
			 G <- symmetrize.graph(g)
			 a <- 1
			 b <- neighborhood(G,a,as.indices=T,open=T)
			 while(T){
			    a <- unique(neighborhood(G,b,as.indices=T,open=T))
			    b <- unique(neighborhood(G,a,as.indices=T,open=T))
				 if(length(c(a,b))>=n) break
			 }
			 if(length(intersect(a,b))>0) stop("g is not bipartite")
			 parts <- rep(0,n)
			 parts[a] <- 1
			 parts[b] <- 2
		 }
		 if(is.list(parts)){
		    x <- 1:n
			 for(i in 1:length(parts)) x[parts[[i]]] <- i
			 parts <- x
		 }
		 if(nrow(d) != ncol(d)){
			 a <- nrow(d)
			 b <- ncol(d)
			 n <- a+b
		    A <- matrix(0,nrow=n,ncol=n)
			 for(i in 1:a){
			    for(j in 1:b){
				    if(d[i,j]){
					    A[i,a+j] <- d[i,j]
					    A[a+j,i] <- d[i,j]
					 }
				 }
			 }
			 d <- A
			 if(!is.null(label)){
			    if(length(label)<n) label <- 1:n
			 }
			 usearrows <- F
			 x <- c(rep(0,a),rep(1,b))
			 y <- c((1:a)/a,(1:b)/b)
		 }
		 else{
			 if(length(parts) != n)
             stop("parts must be a vector of length the number of vertices")
			 uparts <- sort(unique(parts))
			 mp <- max(uparts)
			 if(any(uparts != 1:mp))
			    stop("parts must be a vector containing values from 1 to max(parts)")
		    A <- matrix(0,nrow=length(parts),ncol=length(parts))
			 x <- NULL
			 y <- NULL
			 a <- rep(0,mp)
			 for(i in 1:mp){
			    a[i] <- sum(parts==i)
				 x <- c(x,rep(i,a[i]))
				 y <- c(y,(1:a[i])/a[i])
			 }
			 b <- c(0,cumsum(a))
			 for(i in 1:mp){
			    for(j in 1:mp){
				    if(i != j){
						 A[(b[i]+1):(b[i+1]),(b[j]+1):b[j+1]] <- d[parts==i,parts==j]
					 }
				 }
			 }
			 d <- A
		 }
	 }
    if (jitter) {
		  x <- jitter(x,factor=as.integer(jitter))
		  y <- jitter(y,factor=as.integer(jitter))
    }
	 if(is.null(label.coord)){
		 xl <- x
		 yl <- y
	 }
	 else if(mode != "clique"){
		 xl <- label.coord[,1]
		 yl <- label.coord[,2]
	 }
    use <- displayisolates | !isolated.vertices(g)
	 if(missing(xlim)){
		 xlim <- c(min(x[use]) - pad, max(x[use]) + pad)
		 if(shadow==T && mode=="clique" && keep.circle==T){
		    xlim <- c(min(xlim[1],z-pad),max(xlim[2],z+pad))
		 }
	 }
	 if(missing(ylim)){
		 ylim = c(min(y[use]) - pad, max(y[use]) + pad) 
		 if(shadow==T && mode=="clique" && keep.circle==T){
		    ylim <- c(min(ylim[1],w-pad),max(ylim[2],w+pad))
		 }
	 }
	 if(plotit){
		 if(newplot){
			 if ((length(x) > 0) & (!all(use == FALSE)) & (any(vertex.cex>0))) 
				  plot(x[use], y[use], xlim = xlim, ylim = ylim,
						type = "p", pch = vertex.pch, xlab = xlab,
						ylab = ylab, col = vertex.col, cex = vertex.cex, 
						...)
			 else plot(0, 0, type = "n", pch = vertex.pch, xlab = xlab,
				  ylab = ylab, col = vertex.col, cex = vertex.cex, 
			  ...)
		 }
		 px0 <- vector()
		 py0 <- vector()
		 px1 <- vector()
		 py1 <- vector()
		 e.lwd <- vector()
		 e.col <- vector()
		 e.type <- vector()
		 for (i in c(1:n)[use]) for (j in c(1:n)[use]) if (i != j & d[i,j]>0){
			  px0 <- c(px0, as.real(x[i]))
			  py0 <- c(py0, as.real(y[i]))
			  px1 <- c(px1, as.real(x[j]))
			  py1 <- c(py1, as.real(y[j]))
			  if(is.matrix(edge.lwd)){
			     e.lwd <- c(e.lwd,edge.lwd[i,j])
			  }
			  else if (edge.lwd > 0) 
					e.lwd <- c(e.lwd, edge.lwd)
			  else e.lwd <- c(e.lwd, 1)
			  if(is.matrix(edge.col)){
			     e.col <- c(e.col,edge.col[i,j])
			  }
			  else if (edge.col > 0) 
					e.col <- c(e.col, edge.col)
			  else e.col <- c(e.col, 1)
			  if(is.matrix(edge.type)){
			     e.type <- c(e.type,edge.type[i,j])
			  }
			  else if (edge.type > 0) 
					e.type <- c(e.type, edge.type)
			  else e.type <- c(e.type, 1)
		 }
		 if (usearrows & (length(px0) > 0)) 
			  switch(arrowtype,
			  arrows(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), length = arrowhead.length, angle = 15, 
					col = e.col, lty = e.type, lwd = e.lwd),
			  midarrows(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), length = arrowhead.length, angle = 15, 
					col = e.col, lty = e.type, lwd = e.lwd),
			  endarrows(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), length = arrowhead.length, 
					col = e.col, hcol=head.col, hlwd=head.lwd,
					lty = e.type, lwd = e.lwd))
		 else if (length(px0) > 0) 
			  segments(as.vector(px0), as.vector(py0), as.vector(px1), 
					as.vector(py1), col = e.col, lty = e.type, 
					lwd = e.lwd)
		 if ((length(label) > 0) & (!all(use == FALSE))) 
                      for(i in 1:length(xl)){
			  text(xl[i], yl[i], label[i], pos = 1, cex = label.cex, 
					col = label.col)
                      }
		 if ((length(x) > 0) & (!all(use == FALSE)) & (any(vertex.cex>0))) 
			  points(x[use], y[use], pch = vertex.pch,col = vertex.col, 
			  cex = vertex.cex) 
		 if(mode=="clique" && keep.circle==T && shadow==T){
			 text(z,w,slabel,pos=1,cex=label.cex[1],col=gray(0.8))
		 }
	 }
    invisible(cbind(x,y))
}