Nothing
R/ChowRuskey.R
In Vennerable: Venn and Euler area-proportional diagrams
#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 2, 2019, 5:49 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.