R/ChowRuskey.R

Defines functions deltasmooth deltagivenouter compute.delta raypoint.area raypoint.to.xy compute.CR dual.region.order make.maxiray makeirs makesrp unify.rays make.setlist.from.AWFE make.setlist cycle.to.path remove.loose.nodes node.to.ray scythegr TDtograph

Documented in compute.CR compute.delta deltagivenouter deltasmooth makeirs make.maxiray make.setlist make.setlist.from.AWFE makesrp node.to.ray scythegr TDtograph

#warning("Entering ChowRuskey")

TDtograph <- function(TD) {
	nodes <- names(TD@nodeList)
	# need to delete nonintersection points recognised by two successive 
	# edges in a face corresponding to the same set
	gr <- new("graphNEL",nodes=nodes,edgemode="directed")
	nodeDataDefaults(gr,"x") <- NA; nodeDataDefaults(gr,"y") <- NA; 

	for (node in nodes) {
		nodeData(gr,node,"x") <- TD@nodeList[[node]][,1]
		nodeData(gr,node,"y") <- TD@nodeList[[node]][,2]
	}
	edgeDataDefaults(gr,"Set") <- NA
	edgeDataDefaults(gr,"Name") <- NA
	edgeDataDefaults(gr,"Face+") <- NA
	edgeDataDefaults(gr,"Face-") <- NA
	for (sname in names(TD@setList)) {
		sedges <- .face.to.faceEdges(drawing=TD,face=sname,type="set")
		for (edgeName in names(sedges)) {
			edge <- sedges[[edgeName]]	
			gr <- addEdge(from=edge@from,to=edge@to,graph=gr)
			edgeData(gr,from=edge@from,to=edge@to,"Set") <- which(sname==names(TD@setList))
			edgeData(gr,from=edge@from,to=edge@to,"Name") <- edgeName
		}
	}
	for (fname in .faceNames(TD)) {
		fedgeNames <-  .faceEdgeNames(drawing=TD,face=fname,unsigned=FALSE)
		fedges <- .face.to.faceEdges(drawing=TD,face=fname)
		for (edgeName in fedgeNames ) {
			edge <- fedges[[edgeName ]]
			faceReverse <-  substr(edgeName,1,1)=="-"
			if (!faceReverse) {
				edgeData(gr,from=edge@from,to=edge@to,"Face+") <- fname
			} else {
				edgeData(gr,from=edge@to,to=edge@from,"Face-") <- fname
			}
		}
	}

	gregions <- list()
	for (fname in .faceNames(TD)) {
		fedges <- .face.to.faceEdges(drawing=TD,face=fname)
		fnodes <- sapply(fedges,function(x)x@to)
		gregion <- subGraph(gr,snodes=fnodes)
		gregions[[fname]] <- gregion
	}


	list(graph=gr,regions=gregions)
}

scythegr <- function(TDgr) {
	faces <- unlist(edgeData(TDgr,attr="Face+"))
	vfaces <- setdiff(faces,"DarkMatter")
	nSets <- unique(nchar(vfaces)); stopifnot(length(nSets)==1)
	dualPath <- sapply(0:(nSets),function(i){paste(paste(rep("1",nSets-i),collapse=""),paste(rep("0",i),collapse=""),sep="")})
	dualPath[length(dualPath)] <- "DarkMatter"
	if (nSets==4) {dualPath[2] <- "1101"} # to agree with Chow Ruskey example
	ed <- edgeData(TDgr)
	faces1 <- sapply(ed,function(x)x$"Face+")
	faces2 <- sapply(ed,function(x)x$"Face-")
	for (ix in 1:(length(dualPath)-1)) {
		face1 <- dualPath[ix]; face2 <- dualPath[ix+1]
		edge1 <- faces1[faces1 == face1]
		edge2 <- faces2[faces2 == face2]
		cutedge <- intersect(names(edge1),names(edge2))	; stopifnot(length(cutedge)==1)
		fromto <- strsplit(cutedge,split="|",fixed=TRUE)[[1]]
		TDgr <- removeEdge(graph=TDgr,from=fromto[1],to=fromto[2])
	}
	TDgr
}


node.to.ray <- function(pnodes,atsort) {
	by.tsort <- sapply(pnodes,function(pnode){
		which(sapply(atsort,function(nodes)pnode%in%nodes))}
	)
	names(by.tsort) <- pnodes
	2* (by.tsort)-1
}

.remove.loose.nodes <- function(G) {
	loose.nodes <- nodes(G)[( sapply(edges(G),length)==0 ) & (sapply(inEdges(G),length)==0)]	
	if (length(loose.nodes)>0) {
		G <- removeNode(loose.nodes,G)
	}
	G
}

.cycle.to.path <- function(face) {
	node.path <- nodes(face)[1]
	repeat {
		next.edge <- edges(face,node.path[length(node.path)])[[1]]
		node.path <- c(node.path,next.edge)
		if (next.edge==nodes(face)[1]) { break }
	}
	node.path
}

make.setlist <- function(region,atsort) {
	# we know for each region the fixed and free  edges around it
	# we convert the node names to ray numbers
	# and use the Set attribute to assign ray points to each set
	# on the free boundary
	rays <- node.to.ray(nodes(region),atsort) #the mapping from the topo sort 

	twok <- 2*length(atsort)

	free.region <- region
	fixed.region <- region
	fixed.edges <- unlist(edgeData(region,attr="fixed"))
	for (edgebar in names(fixed.edges)) {
		fromto <- strsplit(edgebar,split="|",fixed=TRUE)[[1]]
		if (fixed.edges[edgebar]) {
			free.region <- removeEdge(fromto[1],fromto[2],free.region)
		} else {
			fixed.region <- removeEdge(fromto[1],fromto[2],fixed.region)
		}
	}
	free.region <- .remove.loose.nodes(free.region) 
	fixed.region <-.remove.loose.nodes(fixed.region) 

	isStart <- (length( fixed.edges[fixed.edges])==0)
	if (isStart) { 
		free.node.path <- .cycle.to.path(free.region)
	} else  {
		free.node.path <- tsort(free.region)
		fixed.node.path <- tsort(fixed.region) 
	}

	EdgeLabels <- lapply(1: (length(free.node.path)-1), function(n) {
		from <- free.node.path[n]
		to <- free.node.path[n+1]
		rayfrom <- rays[from]
		rayto <- rays[to]
		if (rayto<rayfrom) {
			rayrange <- c( (rayfrom:(twok)),1:rayto) 
		} else {
			rayrange <- rayfrom:rayto
		}
		Set <-edgeData(free.region,from=from,to=to,attr="Set")[[1]]
		rayptnames <- rep(NA,length(rayrange))
		list(SetNumber=Set,Raypoints=rayrange,RayPointNames=rays[c(to,from)])	
	})

	if (!isStart) {
		EdgeLabels[[1]]$Raypoints <- EdgeLabels[[1]]$Raypoints[-1] # the first point is on the already drawn set
		EL <- EdgeLabels[[length(EdgeLabels)]]
		EL$Raypoints <- EL$Raypoints[-length(EL$Raypoints) ] # and the last .
		EdgeLabels[[length(EdgeLabels)]] <- EL
	}

	EdgeLabels
}

make.setlist.from.AWFE <- function(G,regions,regionOrder,atsort) {
	setlist <- list()
	edgeDataDefaults(G,"fixed") <- FALSE

	for (vs in regionOrder[-length(regionOrder)]){ # last entry is "0000..."
#cat(vs,"\n");
		region <- regions[[vs]]
		edgeDataDefaults(region,"fixed") <- FALSE
		for (from in names(edges(region))) {
			for (to in edges(region)[[from]]) {
				edgeData(region,from,to,"fixed") <- edgeData(G,from,to,"fixed")
			}
		}
		setlist[[vs]] <- make.setlist(region,atsort)
		
		vs.edges <- edges(region)
		for (from in names(vs.edges)) {
			for (to in edges(region)[[from]]) {
				edgeData(G,from,to,"fixed") <- TRUE
			}
		}
	}
	setlist
}



.unify.rays <- function(asetlines,twok) {
	# given the start and end points of each Set in the rayPoints
	# our task is to join them up and supply the outer two rays
	# (unless we already have the whole circle)
	
	raystarts <- sapply(asetlines,function(x)x$Raypoints[1])
	rayends <- sapply(asetlines,function(x)x$Raypoints[length(x$Raypoints)])
	thru <- rayends < raystarts
	raystartsabs <- ifelse(thru,raystarts-twok,raystarts)
	ors <- raystarts[order(raystartsabs)]
	ore <- rayends[order(raystartsabs)]
	# should all be overlapping 
	stopifnot( length(ors==1) | (ors[2:length(ors)] == ore[1:(length(ore)-1)]))
	rs <- ors[1]
	re <- ore[length(ore)]
	mod2k1 <- function(n) { 1+ (n-1)%% twok }
	if (rs==re) {
		allrays <- 1:twok
	} else {
		rs <- mod2k1(rs-1); 
		re <- mod2k1(re+1);
	
		if (rs>=re) {
			allrays <- c(rs:twok,1:re)
		} else {
			allrays <- (rs:re)
		}
	} 
	allrays
}


makesrp <- function(vs,asetlines,SetRayPoints,Weight,outerRay,angleray) {
	# helper function for compute.CR
	# given an intersection set (vs), with a topology specified in asetlines,
	# first of all compute the (outer) free edges of the intersection set
# { cat(sprintf("%s \n",vs))}
	
	# the rays on which the set is drawn
	raypts <- .unify.rays(asetlines,ncol(SetRayPoints))

	# the offset from the fixed (inner) edges to the free ones
	delta <- deltagivenouter (outerRay[raypts],Weight[vs],angleray)
	# then assign the correct labels to the outer edges
	deltaix <- 1
	for (EL in asetlines) {
		rp <- EL$Raypoints
		# special case for the central point
		if (length(raypts)==ncol(SetRayPoints)) {
			addto <- unique(delta) ; 
			stopifnot(length(addto)==1)
		} else {
			addto <- delta[ seq(deltaix,deltaix+length(rp)-1)]
		}
		SetRayPoints[EL$SetNumber,rp] <- 	outerRay[rp]+addto
		deltaix <- deltaix+length(rp)-1
	}
	SetRayPoints
}

makeirs <- function(vs,asetlines,SetRayPoints,IntersectionRaySets,outerRay) {
	# once the SetRayPoints for the current intersection set have been computed,
	# use them, plus the topology, to compute the perimeter of the intersection set by
	# taking the right bits from the right set

	# the intersection polygon goes out along the inside
	raypts <- .unify.rays(asetlines,ncol(SetRayPoints))
	isStart = (length(raypts)==ncol(SetRayPoints))
	if (!isStart) {
		IntersectionRaySets[[vs]] <- cbind(raypts,outerRay[raypts])
		# then back (hence the rev) along the outside
		for (setline in rev(asetlines)) {
			IntersectionRaySets[[vs]] <- rbind(IntersectionRaySets[[vs]],
				cbind(rev(setline$Raypoints),SetRayPoints[setline$SetNumber,rev(setline$Raypoints)])
			)
		}
	} else {
		IntersectionRaySets[[vs]] <- cbind(numeric(0),numeric(0)) # inner polygon has no inner
		for (setline in (asetlines)) {
			IntersectionRaySets[[vs]] <- rbind(IntersectionRaySets[[vs]],
				cbind((setline$Raypoints),SetRayPoints[setline$SetNumber,(setline$Raypoints)])
			)
		}
	}
	IntersectionRaySets
}

make.maxiray <- function(irs,twok) {
	irs.start <- irs[1,1]
	mod2k1 <- function(n) { 1+ (n-1)%% twok }
	irs[,1] <- mod2k1(irs[,1]-irs[1,1])
	irsrange <- lapply(split(irs[,2],irs[,1]),function(x)(range(x)))
	irsdiff <- sapply(irsrange,function(x)diff((x)))
	if (length(irsdiff)==twok) { # the inner ring
		midray <- 1
		midpoint <- 0
	} else  {
		maxirays <- names(irsdiff)[irsdiff==max(irsdiff)]
		midray <- as.vector(quantile(as.numeric(maxirays),0.5,type=1))
		midpoint <- mean(irsrange[[as.character(midray)]])
		midray <- mod2k1(midray+irs.start)
	}
	res <- matrix(c(midray,midpoint),ncol=2)
	return(res)
}	

.dual.region.order <- function(regionNames ) {
	nnsplit <- strsplit(regionNames ,split="")
	nncount <- sapply(nnsplit,function(x)length(x[x=="1"]))
	regionNames [order(-nncount)]
}

compute.CR <- function(V,doWeights=TRUE) {
	nSets <- NumberOfSets(V)
	Vorig <- V
	if (!doWeights) {	Weights(V) <- 1 + 0 * Weights(V)}
	
	AWFE.diagram <- compute.AWFE(V,type="battle")
	AWFE.graphregions <- TDtograph(AWFE.diagram)
	AWFE.graph <- AWFE.graphregions$graph
	AWFE.regions <- AWFE.graphregions$regions

	sTDgr <- scythegr(AWFE.graph)
	ray.to.point.data <- my.tsort(sTDgr)
	k <- 	length(ray.to.point.data)
	twok <- 2*k

	regionNames <-	.faceNames(AWFE.diagram)
	regionNames[regionNames=="DarkMatter"] <- dark.matter.signature(V)
	regionOrder <- .dual.region.order(regionNames )
	
	####################
	# then prepare the graph we want to construct
	#################
	Weight <- Weights(V)
		
	
	angleray <- 2*pi / twok 
	outerRay <- rep(0,twok )
	SetRayPoints <- matrix(NA,nrow=nSets,ncol=twok )

	# we work in ray points with coordinates measured out along each ray
	# at each stage, the farthest known point out on each ray is in outerRay
	# SetRayPoints, the point at which each Set is on the ray will be built up as we work through 
	# the regions in topological order
	# IntersectionRaySets, is the (Set-labelled) list of ray points specifying the edges of each face
	
	setlines <- make.setlist.from.AWFE (G=AWFE.graph,regions=AWFE.regions,regionOrder,atsort=ray.to.point.data)

	
	# this is the loop in topological order of G*
	for (vs in names(setlines)) { 
#cat(outerRay,"\n")
		SetRayPoints <-  makesrp(vs,setlines[[vs]],SetRayPoints,Weight,outerRay,angleray)		
		outerRay <- apply(SetRayPoints,2,max,na.rm=TRUE)
	}

	srp.df <- melt(SetRayPoints); colnames(srp.df) <- c("SetNumber","Ray","r")
	srp.df$Name <- NA
	# that's all the work done, just encode this as a tissue diagram now
	for (vs in names(setlines)) {
		for (ix in 1:length(setlines[[vs]])) {
			asetlines <- setlines[[vs]][[ix]]
			Set <- asetlines$SetNumber
			rpn <- asetlines$RayPointNames
			for (rp in rpn) {
				srp.df[srp.df$SetNumber==Set & srp.df$Ray ==rp,"Name"] <- names(rpn)[rpn==rp]
			}
		}
	}
	xdf <- .raypoint.to.xy(srp.df[,c("Ray","r")],angleray); colnames(xdf) <- c("x","y")
	srp.df <- cbind(srp.df,xdf)
	points.df <- unique(srp.df[!is.na(srp.df$Name),c("Ray","r","Name","x","y")])	


	# we already have the topology and edge patterns correctly set up in AWFE.diagram
	# all we have to do is adjust the position of its nodes, and redraw its edges
	TD <- AWFE.diagram
	TD@nodeList <- lapply(split(points.df[,c("x","y","Name")],points.df$Name),function(xydf){ w<- data.matrix(xydf[,c("x","y")]);rownames(w)<-xydf$Name;w})

	for (six in 1:length(names(TD@setList))) {
		setName <- names(TD@setList)[[six]]
		sedges <- .face.to.faceEdges(drawing=TD,setName,type="set")
		setray <- srp.df[srp.df$SetNumber==six,]
		setray <- rbind(setray,setray)
		for (sex in 1:length(sedges)) {
			ename <- names(sedges)[sex]
			edge <- sedges[[sex]]
			fromix <- which(setray$Name==edge@from)[1]
			tempray <- setray; tempray$Name[1:fromix]<- NA
			toix <- which(tempray$Name==edge@to)[1]
			xymat <- data.matrix(setray[fromix:toix,c("x","y")])
			edge@xy <- xymat 
			edge@bb <- rbind(apply(xymat,2,min),apply(xymat,2,max))
			# this assumes the previous AWFE edge was a line, which it wouldnt be if we used classic rather than battle AWFE
			TD@edgeList[[ename]] <- edge
		}
	}
	# also delete the empty faces
	emptyFaces <- names(Weight)[Weight==0]
	for (fname in emptyFaces) {
		TD@faceList[[fname]] <- NULL
		TD@faceSignature[[fname]] <- NULL
	}

	TD <- as(TD,"TissueDrawing")
	VD <- new("VennDrawing",TD,Vorig)
	SetLabels <- .default.SetLabelPositions(VD)
	VD <- VennSetSetLabels(VD,SetLabels)
	VD <- .square.universe(VD ,doWeights=doWeights)
	FaceLabels <- .default.FaceLabelPositions(VD)
	VD <- VennSetFaceLabels(VD,FaceLabels)
	VD
}

.raypoint.to.xy <- function(raypoints,angleray) {
	angles <- -angleray * (raypoints[,1])
	xy <- matrix( c(cos(angles), sin(angles)),ncol=2,byrow=FALSE)
	xy <- raypoints[,2] * xy
	xy
}	

.raypoint.area <- function(raypoints,angleray) {
	xy <- .raypoint.to.xy(raypoints,angleray)
	.polygon.area(xy)
}

compute.delta <- function( outervals, wght,angleray) {
	nrays <- length(outervals)
	#ca <- 2
	# try change
	ca <- (nrays-3) 
	singix <- c(1,2,nrays -1,nrays )
	cb <- sum(outervals[singix]) + 2* sum(outervals[-singix])
	cc <- -wght/(sin(angleray)/2)
	delta <- (-cb+sqrt(cb^2-4*ca*cc))/(2*ca)
	delta <- rep(delta,nrays -2)
	delta
}

deltagivenouter <- function( outervals, wght,angleray,smidge=1e-3) {
	# given n outervals as ray points on n successive rays separated by angleray
	# return n-2 deltas corresponding to the middle n-2 rays with area wght
	nray <- length(outervals); stopifnot( nray>2)
	if (wght==0) {
		delta <- rep(0,nray-2)
		return(delta)
	} 
	if (length(unique(outervals[2:(nray-1)]))==1) {
		# a nice uniform inside
		# special case: first time through all the outervals are zero and we have a polygon
		if (all(outervals==0)) {
			delta <- 	sqrt(wght/(length(outervals)*sin(angleray)/2))
			delta <- rep(delta,nray-2)
		} else {
			delta <- compute.delta ( outervals, wght,angleray) 
		}
		return(delta)
	} 
	# not uniform. Can we make it so?
	# where.max <- which(outervals==max(outervals))
	rmax <- max(outervals[-c(1,nray)])
	smidge <- min(smidge,wght/(rmax^2*nray*angleray) )
	rmax <- rmax*(1+smidge)

	feasible.shape.raypoints <- rbind(
		cbind(1:nray,outervals),
		cbind((nray-1):2,rmax))
	feasible.shape.area <- .raypoint.area(feasible.shape.raypoints,angleray )
	if (feasible.shape.area <= wght  ) {
		# that means we can fill in the feasible shape; create a uniform outside and call again to make the outer
		# shape (which will be again uniform ) 
		# nice and smooth; fill in the inner bit
			newouter <- c(outervals[1],rep(rmax,nray-2),outervals[nray])
			delta <- deltagivenouter(newouter, wght- feasible.shape.area,angleray)
			delta <- delta+rmax
			delta <- delta-outervals[2:(nray-1)]
	} else {
		if ( wght < feasible.shape.area  & feasible.shape.area < wght + 2*pi*rmax^2*smidge ) {
			# it's possible that the smidge is too large relative to the weight
			# to allow any feasible solution 
			# if this happens should choose a smaller smidge but we bodge it
			# the effect 
			warning("Try a different smidge; bodging")				
		}
		if (nray>5) {
			# if the maxes occur at the edges, take those out
			if (outervals[2]==max(outervals[-c(1,nray)]))	{
				delta.2 <- smidge
				smidge.triangle <- cbind(c(1,2,2),c(outervals[1],outervals[2]+smidge,outervals[2]))
				smidge.area <- .raypoint.area(smidge.triangle,angleray )
				sub.outer <- outervals[-1]
				sub.outer[1] <- sub.outer[1]+smidge
				sub.delta <- deltagivenouter(sub.outer, wght-smidge.area,angleray)
				delta <- c(delta.2,sub.delta)
				return(delta)
			} else if (outervals[nray-1]==max(outervals[-c(1,nray)])) {
				delta.2 <- smidge
				smidge.triangle <- cbind(c(nray,nray-1,nray-1),c(outervals[nray],outervals[nray-1],outervals[nray-1]+smidge))
				smidge.area <- .raypoint.area(smidge.triangle,angleray )
				sub.outer <- outervals[-nray]
				sub.outer[length(sub.outer)] <- sub.outer[length(sub.outer)]+smidge
				sub.delta <- deltagivenouter(sub.outer, wght-smidge.area,angleray)
				delta <- c(sub.delta,delta.2)
				return(delta)
			}
			# if the max is in the middle, want to divide the weights
			# but lets not bother for now
		} # nray>5

		# mostly though, will need to notch.
		# same delta for each ray (could do much better than this...)
		delta <- compute.delta ( outervals, wght,angleray) 
	}
	return(delta)
}


deltasmooth <- function( outervals, wght,angleray) {
	delta <- compute.delta (outervals, wght,angleray)
	if (length(outervals)!=5) {return(delta)}
	s <- (outervals[1]+outervals[2])/outervals[3]
	if (s>1) {
		delta <- delta * c(1/s,s,1/s)
	}
}

Try the Vennerable package in your browser

Any scripts or data that you put into this service are public.

Vennerable documentation built on May 31, 2017, 4:02 a.m.