R/ICORErandomized.R

In CORE: Cores of Recurrent Events

ICORErandomized <-
function(procno=1,COREobj,boundaries,nprocs=1,rngoffset=0){
minscore<-max(COREobj$minscore,min(COREobj$coreTable[,"score"]))
set.seed(COREobj$seedme)
COREobj$nshuffle<-COREobj$nshuffle-rngoffset
myshuffles<-COREobj$nshuffle%/%nprocs
shuffleres<-COREobj$nshuffle%%nprocs
shuffleskip<-myshuffles*(procno-1)+min(shuffleres,procno-1)
if(procno<=shuffleres)myshuffles<-myshuffles+1
weightList<-vector(mode="list",length=nrow(COREobj$coreTable))
weight<-COREobj$input[,"weight"]
for(i in 1:nrow(COREobj$coreTable)){
	za<-COREobj$input[COREobj$input[,"chrom"]==COREobj$coreTable[i,"chrom"],,drop=F]
	cigt<-za[,"start"]<=COREobj$coreTable[i,"start"]&
		za[,"end"]>=COREobj$coreTable[i,"end"]
	weight[COREobj$input[,"chrom"]==COREobj$coreTable[i,"chrom"]][cigt]<-      
		weight[COREobj$input[,"chrom"]==COREobj$coreTable[i,"chrom"]][cigt]*      
		(1-((COREobj$coreTable[i,"end"]-COREobj$coreTable[i,"start"]+1)/
		(za[cigt,"end"]-za[cigt,"start"]+1))^COREobj$pow)
	weightList[[i]]<-matrix(ncol=2,
		data=c(which(COREobj$input[,"chrom"]==COREobj$coreTable[i,"chrom"])[cigt],
		weight[COREobj$input[,"chrom"]==COREobj$coreTable[i,"chrom"]][cigt]))
}
simscores<-matrix(ncol=myshuffles,nrow=nrow(COREobj$coreTable))
chrmax<-nrow(boundaries)
advanceRNG(randopt=COREobj$shufflemethod,nrand=shuffleskip+rngoffset,
	nevents=nrow(COREobj$input))
for(shuffle in 1:myshuffles){
if(COREobj$shufflemethod=="SIMPLE")z<-
	cbind(randomEventMoves(COREobj$input[,"end"]-COREobj$input[,"start"]+1,
	boundaries),COREobj$input[,"weight"])
if(COREobj$shufflemethod=="RESCALE")z<-
	cbind(randomRescaledEventMoves(COREobj$input[,c("start","end","chrom",
	"weight")],boundaries),COREobj$input[,"weight"])
dimnames(z)[[2]]<-c("start","end","chrom","weight")
chu<-unique(z[,"chrom"])
chc<-rep(0,length(chu))
for(i in 1:length(chu))chc[i]<-sum(z[z[, "chrom"]==chu[i],"weight"])
scoremat<-matrix(nrow=nrow(COREobj$coreTable),ncol=length(chu),data=0)
for(ich in 1:length(chu)){
	wza<-which(z[,"chrom"]==chu[rev(order(chc))][ich])
	za<-z[wza,c("start","end","weight"),drop=F]
	oza<-order(za[,"start"])
	ooza<-order(oza)
	za<-za[oza,,drop=F]
	spairs<-makepairs(za)
	spairs<-spairs[order(spairs[,"start"]),,drop=F]
	psched<-interval.schedule(spairs)
	invl<-1/(za[,"end"]-za[,"start"]+1)
	score<-rep(0,nrow(spairs))
	for(sch in 1:max(psched)){
		pspairs<-spairs[psched==sch,,drop=F]
		ci<-containment.indicator(pspairs[,"start"],pspairs[,"end"],
			za[,"start"],za[,"end"])
		cigt<-(ci[,2]>=ci[,1])
		if(sum(cigt)>0){
			pci<-ci[cigt,,drop=F]
			pinvl<-invl[cigt]
			pwa<-za[cigt,"weight"]
			score[psched==sch]<-score[psched==sch]+
				(spairs[psched==sch,"end"]-
				spairs[psched==sch,"start"]+1)^COREobj$pow*
				fill.values(pci,pwa*pinvl^COREobj$pow,nrow(pspairs))
		}
	}
	for(i in 1:nrow(COREobj$coreTable)){
		spairs<-spairs[score>=minscore,,drop=F]
		if(nrow(spairs)==0){
			scoremat[i:nrow(COREobj$coreTable),ich]<-0
			break
		}
		scoremat[i,ich]<-max(score)
		if(i==nrow(COREobj$coreTable))break
		score<-score[score>=minscore]
		psched<-interval.schedule(spairs)
		whereIwas<-which(wza%in%weightList[[i]][,1])
		if(length(whereIwas)>0){
			whereIamNow<-ooza[whereIwas]
			zaf<-cbind(za[whereIamNow,,drop=F],
				weightList[[i]][weightList[[i]][,1]%in%wza,2])
			dimnames(zaf)[[2]][ncol(zaf)]<-"neweight"
			za[whereIamNow,"weight"]<-zaf[,"neweight"]
			zaf<-zaf[order(zaf[,"start"]),,drop=F]
			invl<-1/(zaf[,"end"]-zaf[,"start"]+1)
			for(sch in 1:max(psched)){
				pspairs<-spairs[psched==sch,,drop=F]
				ci<-containment.indicator(pspairs[,"start"],pspairs[,"end"],
					zaf[,"start"],zaf[,"end"])
				cigt<-(ci[,2]>=ci[,1])
				if(sum(cigt)>0){
					pci<-ci[cigt,,drop=F]
					pinvl<-invl[cigt]
					pwa<-zaf[cigt,"neweight"]-zaf[cigt,"weight"]
					score[psched==sch]<-score[psched==sch]+
						(spairs[psched==sch,"end"]-
						spairs[psched==sch,"start"]+1)^COREobj$pow*
						fill.values(pci,pwa*pinvl^COREobj$pow,nrow(pspairs))
				}
			}
		}
	}
}
simscores[,shuffle]<-apply(scoremat,1,max)
}
return(simscores)
}