##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
getProtocolData.Affy <- function(filenames){
	scanDates <- data.frame(ScanDate=sapply(filenames, celfileDate))
	rownames(scanDates) <- basename(rownames(scanDates))
	protocoldata <- new("AnnotatedDataFrame",
			    data=scanDates,
			    varMetadata=data.frame(labelDescription=colnames(scanDates),
			                           row.names=colnames(scanDates)))
	return(protocoldata)
}
getFeatureData.Affy <- function(cdfName, copynumber=FALSE){
	pkgname <- getCrlmmAnnotationName(cdfName)
	if(!require(pkgname, character.only=TRUE)){
		suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")
		msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)
		message(strwrap(msg))
		stop("Package ", pkgname, " could not be found.")
		rm(suggCall, msg)
	}
	loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)
	loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)
	loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)
	gns <- getVarInEnv("gns")
	path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))
	load(file.path(path, "snpProbes.rda"))

	if(copynumber){
		load(file.path(path, "cnProbes.rda"))
		cnProbes <- get("cnProbes")
		snpIndex <- seq(along=gns)
		npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex) 
		featurenames <- c(gns, rownames(cnProbes))
	} else featurenames <- gns
	fvarlabels=c("chromosome", "position", "isSnp")
	M <- matrix(NA, length(featurenames), 3, dimnames=list(featurenames, fvarlabels))
	index <- match(rownames(snpProbes), rownames(M)) #only snp probes in M get assigned position
	M[index, "position"] <- snpProbes[, grep("pos", colnames(snpProbes))]
	M[index, "chromosome"] <- snpProbes[, grep("chr", colnames(snpProbes))]
	M[index, "isSnp"] <- 1L
	index <- which(is.na(M[, "isSnp"]))
	M[index, "isSnp"] <- 1L

	if(copynumber){
		index <- match(rownames(cnProbes), rownames(M)) #only snp probes in M get assigned position
		M[index, "position"] <- cnProbes[, grep("pos", colnames(cnProbes))]
		M[index, "chromosome"] <- cnProbes[, grep("chr", colnames(cnProbes))]
		M[index, "isSnp"] <- 0L
	}
	##A few of the snpProbes do not match -- I think it is chromosome Y.
	M[is.na(M[, "isSnp"]), "isSnp"] <- 1L
	return(new("AnnotatedDataFrame", data=data.frame(M)))
	##list(snpIndex, npIndex, fns)
	##crlmmOpts$snpRange <- range(snpIndex)
	##crlmmOpts$npRange <- range(npIndex)
}
construct <- function(filenames, cdfName, copynumber=FALSE, sns, verbose=TRUE){
	if(verbose) message("Initializing container for assay data elements alleleA, alleleB, call, callProbability")
	if(missing(sns)) sns <- basename(filenames)
	protocolData <- getProtocolData.Affy(filenames)
	featureData <- getFeatureData.Affy(cdfName, copynumber=copynumber)
	nr <- nrow(featureData); nc <- length(filenames)
	callSet <- new("SnpSuperSet", 
			 alleleA=initializeBigMatrix(name="A", nr, nc),
			 alleleB=initializeBigMatrix(name="B", nr, nc),
			 call=initializeBigMatrix(name="call", nr, nc),
			 callProbability=initializeBigMatrix(name="callPr", nr,nc),
			 featureData=featureData,
			 annotation=cdfName)
	protocolData(callSet) <- protocolData
	sampleNames(callSet) <- sns
	return(callSet)
}
setMethod("close", "AlleleSet", function(con, ...){
	object <- con
	names <- ls(assayData(object))
	L <- length(names)
	for(i in 1:L) close(eval(substitute(assayData(object)[[NAME]], list(NAME=names[i]))))
	return()
})
setMethod("open", "AlleleSet", function(con, ...){
	object <- con
	names <- ls(assayData(object))
	L <- length(names)
	for(i in 1:L) open(eval(substitute(assayData(object)[[NAME]], list(NAME=names[i]))))
	return()
})

##setReplaceMethod("calls", "SnpSuperSet", function(object, value) assayDataElementReplace(object, "call", value))
##setReplaceMethod("confs", "SnpSuperSet", function(object, value) assayDataElementReplace(object, "callProbability", value))
##setMethod("confs", "SnpSuperSet", function(object) assayDataElement(object, "callProbability"))
genotype <- function(filenames, cdfName, mixtureSampleSize=10^5,
		     fitMixture=TRUE,
		     eps=0.1, verbose=TRUE,
		     seed=1, sns, copynumber=FALSE,
		     probs=rep(1/3, 3),
		     DF=6,
		     SNRMin=5,
		     recallMin=10,
		     recallRegMin=1000,
		     gender=NULL,
		     returnParams=TRUE,
		     badSNP=0.7){
	if(missing(cdfName)) stop("must specify cdfName")
	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")
	if(missing(sns)) sns <- basename(filenames)
	## callSet contains potentially very big matrices
	## More big matrices are created within snprma, that will then be removed.
	callSet <- construct(filenames=filenames,
			     cdfName=cdfName,
			     copynumber=copynumber,
			     sns=sns,
			     verbose=verbose)
	mixtureParams <- matrix(NA, 4, length(filenames))
	snp.index <- which(isSnp(callSet)==1)
	batches <- splitIndicesByLength(1:ncol(callSet), ocSamples())
	for(j in batches){
		snprmaRes <- snprma(filenames=filenames[j],
				    mixtureSampleSize=mixtureSampleSize,
				    fitMixture=fitMixture,
				    eps=eps,
				    verbose=verbose,
				    seed=seed,
				    cdfName=cdfName,
				    sns=sns[j])
		stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))
		pData(callSet)$SKW[j] <- snprmaRes$SKW
		pData(callSet)$SNR[j] <- snprmaRes$SNR
		A(callSet)[snp.index, j] <- snprmaRes[["A"]]
		B(callSet)[snp.index, j] <- snprmaRes[["B"]]
		mixtureParams[, j] <- snprmaRes$mixtureParams
		rm(snprmaRes); gc()
	}
	## remove snprmaRes and garbage collect before quantile-normalizing NP probes 
	if(copynumber){
		np.index <- which(isSnp(callSet) == 0)
		cnrmaRes <- cnrma(filenames=filenames[j],
				  cdfName=cdfName,
				  row.names=featureNames(callSet)[np.index],				  
				  sns=sns[j],
				  seed=seed,
				  verbose=verbose)
		stopifnot(identical(featureNames(callSet)[np.index], rownames(cnrmaRes)))
		A(callSet)[np.index, j] <- cnrmaRes
		rm(cnrmaRes); gc()
	}
	for(j in batches){
		tmp <- crlmmGT(A=A(callSet)[snp.index, j],
			       B=B(callSet)[snp.index, j],
			       SNR=callSet$SNR[j],
			       mixtureParams=mixtureParams[, j],
			       cdfName=annotation(callSet),
			       row.names=featureNames(callSet)[snp.index],
			       col.names=sampleNames(callSet)[j],
			       probs=probs,
			       DF=DF,
			       SNRMin=SNRMin,
			       recallMin=recallMin,
			       recallRegMin=recallRegMin,
			       gender=gender,
			       verbose=verbose,
			       returnParams=returnParams,
			       badSNP=badSNP)
		snpCall(callSet)[snp.index, j] <- tmp[["calls"]]
		snpCallProbability(callSet)[snp.index, j] <- tmp[["confs"]]
		callSet$gender[j] <- tmp$gender
		rm(tmp); gc()
	}	
	return(callSet)
}

##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
## For Illumina
##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
getPhenoData <- function(sampleSheet=NULL, arrayNames=NULL,
			 arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A")){
	if(!is.null(arrayNames)) {
		pd = new("AnnotatedDataFrame", data = data.frame(Sample_ID=arrayNames))
	}
	if(!is.null(sampleSheet)) { # get array info from Illumina's sample sheet
		if(is.null(arrayNames)){
			##arrayNames=NULL
			if(!is.null(arrayInfoColNames$barcode) && (arrayInfoColNames$barcode %in% colnames(sampleSheet))) {
				barcode = sampleSheet[,arrayInfoColNames$barcode]
				arrayNames=barcode
			}
			if(!is.null(arrayInfoColNames$position) && (arrayInfoColNames$position %in% colnames(sampleSheet))) {  
				position = sampleSheet[,arrayInfoColNames$position]
				if(is.null(arrayNames))
					arrayNames=position
				else
					arrayNames = paste(arrayNames, position, sep=sep)
				if(highDensity) {
					hdExt = list(A="R01C01", B="R01C02", C="R02C01", D="R02C02")
					for(i in names(hdExt))
						arrayNames = sub(paste(sep, i, sep=""), paste(sep, hdExt[[i]], sep=""), arrayNames)
				}
			}
		}
		pd = new("AnnotatedDataFrame", data = sampleSheet)
		sampleNames(pd) <- basename(arrayNames)               
	}
	if(is.null(arrayNames)) {
		arrayNames = gsub(paste(sep, fileExt$green, sep=""), "", dir(pattern=fileExt$green, path=path))
		if(!is.null(sampleSheet)) {
			sampleSheet=NULL
			cat("Could not find required info in \'sampleSheet\' - ignoring.  Check \'sampleSheet\' and/or \'arrayInfoColNames\'\n")
		}
		pd = new("AnnotatedDataFrame", data = data.frame(Sample_ID=arrayNames))
	}
	return(pd)
}
constructRG <- function(filenames, cdfName, sns, verbose, fileExt, sep, sampleSheet, arrayInfoColNames){
	if(verbose)	message("reading first idat file to extract feature data")
	grnfile <- paste(filenames[1], fileExt$green, sep=sep)
	if(!file.exists(grnfile)){
                stop(paste(grnfile, " does not exist. Check fileExt argument"))
        }
        G <- readIDAT(grnfile)
        idsG = rownames(G$Quants)
        nr <- length(idsG)
	fD <- new("AnnotatedDataFrame", data=data.frame(row.names=idsG))##, varMetadata=data.frame(labelDescript
	nr <- nrow(fD)
	dns <- list(featureNames(fD), basename(filenames))
	RG <- new("NChannelSet",
		  R=initializeBigMatrix(name="R", nr=nr, nc=length(filenames)),
		  G=initializeBigMatrix(name="G", nr=nr, nc=length(filenames)),
		  zero=initializeBigMatrix(name="zero", nr=nr, nc=length(filenames)),
		  featureData=fD,
		  annotation=cdfName)
	phenoData(RG) <- getPhenoData(sampleSheet=sampleSheet, arrayNames=filenames,
				      arrayInfoColNames=arrayInfoColNames)
##	pD <- data.frame(matrix(NA, length(sampleNames(RG)), 12), row.names=sampleNames(RG))
##	colnames(pD) <- c("Index","HapMap.Name","Name","ID",
##			  "Gender", "Plate", "Well", "Group", "Parent1",
##			  "Parent2","Replicate","SentrixPosition")
	##phenoData(RG) <- new("AnnotatedDataFrame", data=pD)
	pD <- data.frame(matrix(NA, length(sampleNames(RG)), 1), row.names=sampleNames(RG))
	colnames(pD) <- "ScanDate"
	protocolData(RG) <- new("AnnotatedDataFrame", data=pD)
	sampleNames(RG) <- basename(filenames)
	storageMode(RG) <- "environment"
	RG##featureData=ops$illuminaOpts[["featureData"]])
}
crlmmIlluminaRS <- function(sampleSheet=NULL,
			    arrayNames=NULL,
			    ids=NULL,
			    path=".",
			    arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A"),
			    highDensity=FALSE,
			    sep="_",
			    fileExt=list(green="Grn.idat", red="Red.idat"),
			    stripNorm=TRUE,
			    useTarget=TRUE,
			    row.names=TRUE, 
			    col.names=TRUE,
			    probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,
			    seed=1, save.ab=FALSE, snpFile, cnFile,
			    mixtureSampleSize=10^5, eps=0.1, verbose=TRUE,
			    cdfName, sns, recallMin=10, recallRegMin=1000,
			    returnParams=FALSE, badSNP=.7) {
	if(missing(cdfName)) stop("must specify cdfName")
	if(!isValidCdfName(cdfName)) stop("cdfName not valid.  see validCdfNames")
	if(missing(sns)) sns <- basename(arrayNames)
	RG <- constructRG(filenames=arrayNames,
			  cdfName=cdfName,
			  sns=sns,
			  verbose=verbose,
			  fileExt=fileExt,
			  sep=sep,
			  sampleSheet=sampleSheet,
			  arrayInfoColNames=arrayInfoColNames)
	batches <- splitIndicesByLength(seq(along=arrayNames), ocSamples())
	for(j in batches){
		tmp <- readIdatFiles(sampleSheet=sampleSheet[j, ],
				     arrayNames=arrayNames[j],
				     ids=ids,
				     path=path,
				     arrayInfoColNames=arrayInfoColNames,
				     highDensity=highDensity,
				     sep=sep,
				     fileExt=fileExt,
				     saveDate=TRUE)
		if(nrow(tmp) != nrow(RG)){
			##RS: I don't understand why the IDATS for the
			##same platform do not have necessarily have
			##the same length
			tmp <- tmp[featureNames(tmp) %in% featureNames(RG), ]
		}
		index <- match(featureNames(tmp), featureNames(RG))
		assayData(RG)$R[index,j] <- assayData(tmp)$R
		assayData(RG)$G[index, j] <- assayData(tmp)$G
		assayData(RG)$zero[index, j] <- assayData(tmp)$zero
		pData(RG)[j, ] <- pData(tmp)
		annotation(RG) <- annotation(tmp)
		pData(protocolData(RG))[j, ] <- pData(protocolData(tmp))
		rm(tmp); gc()
	}
	message("Saving RG")
	save(RG, file=file.path(ldPath(), "RG.rda"))
	k <- 1
	for(j in batches){
		##MR: Might want to to nonpolymorphic markers in a
		##separate step and make that optional
		tmp <- RGtoXY(RG[, j], chipType=cdfName)
		if(k == 1){
			##initialize XYSet
			nc <- ncol(tmp)
			nr <- nrow(tmp)
			XY <- new("NChannelSet",
				  X=initializeBigMatrix(name="X", nr, nc),
				  Y=initializeBigMatrix(name="Y", nr, nc),
				  zero=initializeBigMatrix(name="zero", nr, nc),
				  experimentData=experimentData(RG),
				  phenoData=phenoData(RG),
				  protocolData=protocolData(RG),
				  annotation=annotation(RG))
		}
		storageMode(XY) <- "environment"
		assayData(XY)$X[, j] <- assayData(tmp)$X
		assayData(XY)$Y[, j] <- assayData(tmp)$Y
		assayData(XY)$zero[, j] <- assayData(tmp)$zero
		k <- k+1
		rm(tmp); gc()
	}
	annotation(XY) <- cdfName
	message("Saving XY")
	save(XY, file=file.path(ldPath(), "XY.rda"))
	mixtureParams <- matrix(NA, 4, length(filenames))
	k <- 1
	for(j in batches){
		res <- preprocessInfinium2(XY[, j],
					   mixtureSampleSize=mixtureSampleSize,
					   fitMixture=TRUE,
					   verbose=verbose,
					   seed=seed,
					   eps=eps,
					   cdfName=cdfName,
					   sns=sns[j],
					   stripNorm=stripNorm,
					   useTarget=useTarget)
		if(k==1){
			## MR: number of rows should be number of SNPs + number of nonpolymorphic markers.
			##  Here, I'm just using the # of rows returned from the above function
			callSet <- new("SnpSuperSet",
				       alleleA=initializeBigMatrix(name="A", nr=nrow(res[[1]]), nc=ncol(XY)),
				       alleleB=initializeBigMatrix(name="B", nr=nrow(res[[2]]), nc=ncol(XY)),
				       call=initializeBigMatrix(name="call", nr=nrow(res[[1]]), nc=ncol(XY)),
				       callProbability=initializeBigMatrix(name="callPr", nr=nrow(res[[1]]), nc=ncol(XY)),
				       protocolData=protocolData(XY),
				       experimentData=experimentData(XY),
				       phenoData=phenoData(XY),
				       annotation=annotation(XY))
			featureNames(callSet) <- res[["gns"]]
			sampleNames(callSet) <- sns
		}
		## MR: we need to define a snp.index
		A(callSet)[, j] <- res[["A"]]
		B(callSet)[, j] <- res[["B"]]
		pData(callSet)$SKW[j] <- res$SKW
		pData(callSet)$SNR[j] <- res$SNR
		mixtureParams[, j] <- res$mixtureParams
		rm(res); gc()
	}
	message("Saving callSe")
	save(callSet, file=file.path(ldPath(), "callSet.rda"))	
	for(j in batches){
		##MR:  edit snp.index
		snp.index <- 1:nrow(callSet)
		tmp <- crlmmGT(A=A(callSet)[snp.index, j],
			       B=B(callSet)[snp.index, j],
			       SNR=callSet$SNR[j],
			       mixtureParams=mixtureParams[, j],
			       cdfName=annotation(callSet),
			       row.names=featureNames(callSet)[snp.index],
			       col.names=sampleNames(callSet)[j],
			       probs=probs,
			       DF=DF,
			       SNRMin=SNRMin,
			       recallMin=recallMin,
			       recallRegMin=recallRegMin,
			       gender=gender,
			       verbose=verbose,
			       returnParams=returnParams,
			       badSNP=badSNP)
		snpCall(callSet)[snp.index, j] <- tmp[["calls"]]
		## many zeros here (?)
		snpCallProbability(callSet)[snp.index, j] <- tmp[["confs"]]
		callSet$gender[j] <- tmp$gender
		rm(tmp); gc()
	}
	return(callSet)
}
##---------------------------------------------------------------------------
##---------------------------------------------------------------------------

rowCovs <- function(x, y, ...){
	notna <- !is.na(x)
	N <- rowSums(notna)
	x <- suppressWarnings(log2(x))
	if(!missing(y)) y <- suppressWarnings(log2(y))
	return(rowSums((x - rowMeans(x, ...)) * (y - rowMeans(y, ...)), ...)/(N-1))
}

rowMAD <- function(x, y, ...){
	notna <- !is.na(x)
	mad <- 1.4826*rowMedians(abs(x-rowMedians(x, ...)), ...)
	return(mad)
}

rowCors <- function(x, y, ...){
	N <- rowSums(!is.na(x))
	x <- suppressWarnings(log2(x))
	y <- suppressWarnings(log2(y))
	sd.x <- rowSds(x, ...)
	sd.y <- rowSds(y, ...)
	covar <- rowSums((x - rowMeans(x, ...)) * (y - rowMeans(y, ...)), ...)/(N-1)
	return(covar/(sd.x*sd.y))
}

generateX <- function(w, X) as.numeric(diag(w) %*% X)
generateIXTX <- function(x, nrow=3) {
	X <- matrix(x, nrow=nrow)
	XTX <- crossprod(X)
	solve(XTX)
}


nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
	I <- isSnp(object)
	Ystar <- Ystar[I, ]
	rownames(Ystar) <- featureNames(object)[isSnp(object)]
	complete <- rowSums(is.na(W)) == 0 & I
	W <- W[I, ]
	if(any(!is.finite(W))){## | any(!is.finite(V))){
		i <- which(rowSums(!is.finite(W)) > 0)
		##browser()
		stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
	}	
	Ns <- tmp.objects[["Ns"]]	
	Ns <- Ns[I, ]
	
	CHR <- unique(chromosome(object))
	batch <- unique(object$batch)
	nuA <- getParam(object, "nuA", batch)
	nuB <- getParam(object, "nuB", batch)
	nuA.se <- getParam(object, "nuA.se", batch)
	nuB.se <- getParam(object, "nuB.se", batch)
	phiA <- getParam(object, "phiA", batch)
	phiB <- getParam(object, "phiB", batch)
	phiA.se <- getParam(object, "phiA.se", batch)
	phiB.se <- getParam(object, "phiB.se", batch)	
	if(CHR==23){
		phiAX <- getParam(object, "phiAX", batch)
		phiBX <- getParam(object, "phiBX", batch)		
	}
	NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05
	if(missing(allele)) stop("must specify allele")
	if(CHR == 23){
		##Design matrix for X chromosome depends on whether there was a sufficient number of AB genotypes
		if(length(grep("AB", colnames(W))) > 0){
			if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
			if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
		} else{
			if(allele == "A") X <- cbind(1, c(1, 0, 2, 0), c(0, 1, 0, 2))
			if(allele == "B") X <- cbind(1, c(0, 1, 0, 2), c(1, 0, 2, 0))
		}
	} else {##autosome
		if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
		if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous		
	}
	##How to quickly generate Xstar, Xstar = diag(W) %*% X
	Xstar <- apply(W, 1, generateX, X)
	IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
	##as.numeric(diag(W[1, ]) %*% X)
	if(CHR == 23){
		betahat <- matrix(NA, 3, nrow(Ystar))
		ses <- matrix(NA, 3, nrow(Ystar))		
	} else{
		betahat <- matrix(NA, 2, nrow(Ystar))
		ses <- matrix(NA, 2, nrow(Ystar))
	}
	for(i in 1:nrow(Ystar)){
		betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
		ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
		ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
	}
	if(allele == "A"){
		nuA[complete] <- betahat[1, ]
		phiA[complete] <- betahat[2, ]
		nuA.se[complete] <- ses[1, ]
		phiA.se[complete] <- ses[2, ]
		object <- pr(object, "nuA", batch, nuA)
		object <- pr(object, "nuA.se", batch, nuA.se)
		object <- pr(object, "phiA", batch, phiA)
		object <- pr(object, "phiA.se", batch, phiA.se)
		if(CHR == 23){
			phiAX[complete] <- betahat[3, ]
			object <- pr(object, "phiAX", batch, phiAX)			
		}
	}
	if(allele=="B"){
		nuB[complete] <- betahat[1, ]
		phiB[complete] <- betahat[2, ]
		nuB.se[complete] <- ses[1, ]
		phiB.se[complete] <- ses[2, ]
		object <- pr(object, "nuB", batch, nuB)
		object <- pr(object, "nuB.se", batch, nuB.se)
		object <- pr(object, "phiB", batch, phiB)		
		object <- pr(object, "phiB.se", batch, phiB.se)
		if(CHR == 23){
			phiBX[complete] <- betahat[3, ]
			object <- pr(object, "phiBX", batch, phiBX)
		}
	}
##	THR.NU.PHI <- cnOptions$THR.NU.PHI
##	if(THR.NU.PHI){
##		verbose <- cnOptions$verbose
##		if(verbose) message("Thresholding nu and phi")
##		object <- thresholdModelParams(object, cnOptions)
##	}
	return(object)
}



predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){
	pkgname <- paste(cdfName, "Crlmm", sep="")
	path <- system.file("extdata", package=pkgname)
        loader("23index.rda", .crlmmPkgEnv, pkgname)
	index <- getVarInEnv("index")
	##load(file.path(path, "23index.rda"))
	XIndex <- index[[1]]	
	if(class(res) != "ABset"){
		A <- res$A
		B <- res$B
	} else{
		A <- res@assayData[["A"]]
		B <- res@assayData[["B"]]
	}
	tmp <- which(rowSums(is.na(A)) > 0)
	XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median, na.rm=TRUE)/2
	SNR <- res$SNR
	if(sum(SNR>SNRMin)==1){
		gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5)))
	}else{
		gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]
	}
	return(gender)
}

##combineIntensities <- function(res, NP=NULL, callSet){
##	rownames(res$B) <- rownames(res$A) <- res$gns
##	colnames(res$B) <- colnames(res$A) <- res$sns
##	if(!is.null(NP)){
##		blank <- matrix(NA, nrow(NP), ncol(NP))
##		dimnames(blank) <- dimnames(NP)
##		A <- rbind(res$A, NP)
##		B <- rbind(res$B, blank)
##	} else {
##		A <- res$A
##		B <- res$B
##	}
##	dimnames(B) <- dimnames(A)
##	index.snps <- match(res$gns, rownames(A))
##	callsConfs <- calls <- matrix(NA, nrow(A), ncol(A), dimnames=dimnames(A))
##	
##	calls[index.snps, ] <- calls(callSet)
##	callsConfs[index.snps, ] <- assayData(callSet)[["callProbability"]]
##	fd <- data.frame(matrix(NA, nrow(calls), length(fvarLabels(callSet))))
##	fd[index.snps, ] <- fData(callSet)
##	rownames(fd) <- rownames(A)
##	colnames(fd) <- fvarLabels(callSet)
##	fD <- new("AnnotatedDataFrame",
##		  data=data.frame(fd),
##		  varMetadata=data.frame(labelDescription=colnames(fd), row.names=colnames(fd)))
##	superSet <- new("CNSet",
##			CA=matrix(NA, nrow(A), ncol(A), dimnames=dimnames(A)),
##			CB=matrix(NA, nrow(A), ncol(A), dimnames=dimnames(A)),
##			alleleA=A,
##			alleleB=B,
##			call=calls,
##			callProbability=callsConfs,
##			featureData=fD,
##			phenoData=phenoData(callSet),
##			experimentData=experimentData(callSet),
##			protocolData=protocolData(callSet),
##			annotation=annotation(callSet))
##	return(superSet)
##}
##
harmonizeDimnamesTo <- function(object1, object2){
	#object2 should be a subset of object 1
	object2 <- object2[featureNames(object2) %in% featureNames(object1), ]
	object1 <- object1[match(featureNames(object2), featureNames(object1)), ]
	object1 <- object1[, match(sampleNames(object2), sampleNames(object1))]
	stopifnot(all.equal(featureNames(object1), featureNames(object2)))
	stopifnot(all.equal(sampleNames(object1), sampleNames(object2)))
	return(object1)
}



crlmmCopynumber <- function(object,
			    batch,
			    chromosome=1:23,
			    MIN.SAMPLES=10,
			    SNRMin=5,
			    MIN.OBS=3,
			    DF.PRIOR=50,
			    bias.adj=FALSE,
			    prior.prob=rep(1/4,4),
			    seed=1,
			    verbose=TRUE,
			    GT.CONF.THR=0.99,
			    PHI.THR=2^6,
			    nHOM.THR=5,
			    MIN.NU=2^3,
			    MIN.PHI=2^3,
			    THR.NU.PHI=TRUE,
			    thresholdCopynumber=TRUE){
	which.batch <- cnOptions[["whichbatch"]]
	cnSet <- new("CNSetLM",
		     alleleA=A(object),
		     alleleB=B(object),
		     call=snpCall(object),
		     callProbability=snpCallProbability(object),
		     CA=initializeBigMatrix("CA", nrow(object), ncol(object)),
		     CB=initializeBigMatrix("CB", nrow(object), ncol(object)),
		     annotation=annotation(object),
		     featureData=featureData(object),
		     experimentData=experimentData(object),
		     phenoData=phenoData(object))
	lM(cnSet) <- initializeParamObject(list(featureNames(cnSet), unique(cnSet$batch)))
	rm(object); gc()
	if(any(cnSet$SNR < SNRMin)){
		message("Excluding ", sum(cnSet$SNR < SNRMin), " samples with SNR below ", SNRMin)
		cnSet <- cnSet[, cnSet$SNR > SNRMin]
	}
	batches <- split(1:ncol(cnSet), cnSet$batch)
	if(any(sapply(batches, length) < MIN.SAMPLES)) message("Excluding batches with fewer than ", MIN.SAMPLES, " samples")
	batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
	for(i in chromosome){
		cat("Chromosome ", i, "\n")
		if(i >= 24) next()
		for(j in batches){
			row.index <- which(chromosome(cnSet) == i)
			tmp <- cnSet[row.index, j]
			featureData(tmp) <- lm.parameters(tmp)
			tmp <- computeCopynumber(tmp,
						 SNRMin=SNRMin,
						 MIN.OBS=MIN.OBS,
						 DF.PRIOR=DF.PRIOR,
						 bias.adj=bias.adj,
						 prior.prob=prior.prob,
						 seed=seed,
						 verbose=verbose,
						 GT.CONF.THR=GT.CONF.THR,
						 PHI.THR=PHI.THR,
						 nHOM.THR=nHOM.THR,
						 MIN.NU=MIN.NU,
						 MIN.PHI=MIN.PHI,
						 thresholdCopynumber=thresholdCopynumber)
			fData(tmp) <- fData(tmp)[, -(1:3)]
			CA(cnSet)[row.index, j] <- tmp@assayData[["CA"]]
			CB(cnSet)[row.index, j] <- tmp@assayData[["CB"]]
			labels.asis <- fvarLabels(tmp)
			labels.asis <- gsub("_", ".", labels.asis)
			k <- match(labels.asis, colnames(lM(cnSet)))
			lM(cnSet)[row.index, k] <- fData(tmp)
			rm(tmp); gc()
		}
	}
	return(cnSet)
}



crlmmWrapper <- function(filenames, cnOptions, ...){
	cdfName <- cnOptions[["cdfName"]]
	load.it <- cnOptions[["load.it"]]
	save.it <- cnOptions[["save.it"]]
	splitByChr <- cnOptions[["splitByChr"]]
	crlmmFile <- cnOptions[["crlmmFile"]]
	intensityFile=cnOptions[["intensityFile"]]
	rgFile=cnOptions[["rgFile"]]
	##use.ff=cnOptions[["use.ff"]]
	outdir <- cnOptions[["outdir"]]
	if(missing(cdfName)) stop("cdfName is missing -- a valid cdfName is required.  See crlmm:::validCdfNames()")
	platform <- whichPlatform(cdfName)
	if(!(platform %in% c("affymetrix", "illumina"))){
		stop("Only 'affymetrix' and 'illumina' platforms are supported at this time.")
	} else {
		message("Checking whether annotation package for the ", platform, " platform is available")
		if(!isValidCdfName(cdfName)){
			cat("FALSE\n")
			stop(cdfName, " is not a valid entry.  See crlmm:::validCdfNames(platform)")
		} else cat("TRUE\n")
	}
	if(missing(intensityFile)) stop("must specify 'intensityFile'.")
	if(missing(crlmmFile)) stop("must specify 'crlmmFile'.")
	if(platform == "illumina"){
		if(missing(rgFile)){
			##stop("must specify 'rgFile'.")
			rgFile <- file.path(dirname(crlmmFile), "rgFile.rda")
			message("rgFile not specified.  Using ", rgFile)
		}
		if(!load.it){
			RG <- readIdatFiles(...)
			if(save.it) save(RG, file=rgFile)
		}
		if(load.it & !file.exists(rgFile)){
			message("load.it is TRUE, but rgFile not present.  Attempting to read the idatFiles.")
			RG <- readIdatFiles(...)
			if(save.it) save(RG, file=rgFile)
		}
		if(load.it & file.exists(rgFile)){
			message("Loading RG file")
			load(rgFile)
			RG <- get("RG")
		}
	}
	if(!(file.exists(dirname(crlmmFile)))) stop(dirname(crlmmFile), " does not exist.")
	if(!(file.exists(dirname(intensityFile)))) stop(dirname(intensityFile), " does not exist.")
	##---------------------------------------------------------------------------
	## FIX
	outfiles <- file.path(dirname(crlmmFile), paste("crlmmSetList_", 1:24, ".rda", sep=""))
	if(load.it & all(file.exists(outfiles))){
		load(outfiles[1])
		crlmmSetList <- get("crlmmSetList")
		if(!all(sampleNames(crlmmSetList) == basename(filenames))){
			stop("load.it is TRUE, but sampleNames(crlmmSetList != basename(filenames))")
		} else{
			return("load.it is TRUE and 'crlmmSetList_<CHR>.rda' objects found. Nothing to do...")
		}
	}
	if(load.it){
		if(!file.exists(crlmmFile)){
			message("load.it is TRUE, but ", crlmmFile, " does not exist.  Rerunning the genotype calling algorithm") 
			load.it <- FALSE
		}
	}
	if(platform == "affymetrix"){
		if(!file.exists(crlmmFile) | !load.it){
			callSet <- crlmm(filenames=filenames,
					     cdfName=cdfName,
					     save.it=TRUE,
					     load.it=load.it,
					     intensityFile=intensityFile)
			message("Quantile normalizing the copy number probes...")		
			cnrmaResult <- cnrma(filenames=filenames, cdfName=cdfName, outdir=outdir)
			if(save.it){
				message("Saving callSet and cnrmaResult to", crlmmFile)
				save(callSet, cnrmaResult, file=crlmmFile)
			}
		} else {
			message("Loading ", crlmmFile, "...")
			load(intensityFile)				
			load(crlmmFile)
			callSet <- get("callSet")
			cnrmaResult <- get("cnrmaResult")
		}
		scanDates <- data.frame(ScanDate=sapply(filenames, celfileDate))
		protocolData(callSet) <- new("AnnotatedDataFrame",
					     data=scanDates,
					     varMetadata=data.frame(labelDescription=colnames(scanDates),
					     row.names=colnames(scanDates)))
	}
	if(platform == "illumina"){
		if(!file.exists(crlmmFile) | !load.it){		
			callSet <- crlmmIllumina(RG=RG,
						 cdfName=cdfName,
						 sns=sampleNames(RG),
						 returnParams=TRUE,
						 save.it=TRUE,
						 intensityFile=intensityFile)
			if(save.it) save(callSet, file=crlmmFile)
		} else {
			message("Loading ", crlmmFile, "...")
			load(crlmmFile)
			callSet <- get("callSet")
		}
		protocolData(callSet) <- protocolData(RG)
	}
	if(platform=="affymetrix") {
		protocolData(callSet)[["ScanDate"]] <- as.character(celDates(filenames))
		sampleNames(protocolData(callSet)) <- sampleNames(callSet)
	}	
	load(intensityFile)
	snprmaResult <- get("res")
	if(platform=="illumina"){
		if(exists("cnAB")){
			np.A <- as.integer(cnAB$A)
			np.B <- as.integer(cnAB$B)
			np <- ifelse(np.A > np.B, np.A, np.B)
			np <- matrix(np, nrow(cnAB$A), ncol(cnAB$A))
			rownames(np) <- cnAB$gns
			colnames(np) <- cnAB$sns
			cnAB$NP <- np
			##sampleNames(callSet) <- res$sns
			sampleNames(callSet) <- cnAB$sns
			cnrmaResult <- get("cnAB")
		} else cnrmaResult <- NULL
	}
	if(platform=="affymetrix"){
		if(exists("cnrmaResult")){
			cnrmaResult <- get("cnrmaResult")
		} else cnrmaResult <- NULL
	}
	crlmmResults <- list(snprmaResult=snprmaResult,
			     cnrmaResult=cnrmaResult,
			     callSet=callSet)

	if(!save.it){
		message("Cleaning up")		
		unlink(intensityFile)
	}
	return(crlmmResults)
}

validCdfNames <- function(){
	c("genomewidesnp6",
	  "genomewidesnp5",
	  "human370v1c",
	  "human370quadv3c",
	  "human550v3b",
	  "human650v3a",
	  "human610quadv1b",
	  "human660quadv1a",
	  "human1mduov3b")
}

isValidCdfName <- function(cdfName){
	chipList <- validCdfNames()
	result <- cdfName %in% chipList	
	if(!(result)){
		warning("cdfName must be one of the following: ",
			chipList)
	}
	return(result)
}

whichPlatform <- function(cdfName){
	index <- grep("genomewidesnp", cdfName)
	if(length(index) > 0){
		platform <- "affymetrix"
	} else{
		index <- grep("human", cdfName)
		platform <- "illumina"
	}
	return(platform)
}


# steps: quantile normalize hapmap: create 1m_reference_cn.rda object
cnrma <- function(filenames, cdfName, row.names, sns, seed=1, verbose=FALSE){
	if(missing(cdfName)) stop("must specify cdfName")
	pkgname <- getCrlmmAnnotationName(cdfName)
	require(pkgname, character.only=TRUE) || stop("Package ", pkgname, " not available")
	if (missing(sns)) sns <- basename(filenames)
        loader("npProbesFid.rda", .crlmmPkgEnv, pkgname)
	fid <- getVarInEnv("npProbesFid")
	fid <- fid[match(row.names, names(fid))]
	set.seed(seed)
	idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything
	SKW <- vector("numeric", length(filenames))
	NP <- matrix(NA, length(fid), length(filenames))
	verbose <- TRUE
	if(verbose){
		message("Processing ", length(filenames), " files.")
		if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)
	}
	if(cdfName=="genomewidesnp6"){
		loader("1m_reference_cn.rda", .crlmmPkgEnv, pkgname)
	}
	if(cdfName=="genomewidesnp5"){
		loader("5.0_reference_cn.rda", .crlmmPkgEnv, pkgname)
	}
	reference <- getVarInEnv("reference")
	##if(!is.matrix(reference)) stop("target distribution for quantile normalization not available.")
	for(i in seq(along=filenames)){
		y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])
		x <- log2(y[idx2])
		SKW[i] <- mean((x-mean(x))^3)/(sd(x)^3)
		rm(x)
		NP[, i] <- as.integer(normalize.quantiles.use.target(y, target=reference))
		if (verbose)
			if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)
			else cat(".")
	}
	dimnames(NP) <- list(names(fid), sns)
	##dimnames(NP) <- list(map[, "man_fsetid"], sns)
	##res3 <- list(NP=NP, SKW=SKW)
	cat("\n")
	return(NP)
}

getFlags <- function(object, PHI.THR){
	batch <- unique(object$batch)
	nuA <- getParam(object, "nuA", batch)
	nuB <- getParam(object, "nuB", batch)
	phiA <- getParam(object, "phiA", batch)
	phiB <- getParam(object, "phiB", batch)
	negativeNus <- nuA < 1 | nuB < 1
	negativePhis <- phiA < PHI.THR | phiB < PHI.THR
	negativeCoef <- negativeNus | negativePhis
	notfinitePhi <- !is.finite(phiA) | !is.finite(phiB)
	flags <- negativeCoef | notfinitePhi
	return(flags)
}


instantiateObjects <- function(object, cnOptions){
	Ns <- matrix(NA, nrow(object), 5)
	colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
	vA <- vB <- muB <- muA <- Ns
	normal <- matrix(TRUE, nrow(object), ncol(object))
	dimnames(normal) <- list(featureNames(object), sampleNames(object))
	tmp.objects <- list(vA=vA,
			    vB=vB,
			    muB=muB,
			    muA=muA,
			    Ns=Ns,
			    normal=normal)
        return(tmp.objects)
}

thresholdCopynumber <- function(object){
	ca <- CA(object)
	cb <- CB(object)
	ca[ca < 0.05] <- 0.05
	ca[ca > 5] <- 5
	cb[cb < 0.05] <- 0.05
	cb[cb > 5] <- 5
	CA(object) <- ca
	CB(object) <- cb
	return(object)
}

cnOptions <- function(
		      MIN.OBS=3,
		      DF.PRIOR=50,
		      bias.adj=FALSE,
		      prior.prob=rep(1/4, 4),
		      seed=123,
		      verbose=TRUE,
		      GT.CONF.THR=0.99,
		      PHI.THR=2^6,##used in nonpolymorphic fxn for training
		      nHOM.THR=5, ##used in nonpolymorphic fxn for training
		      MIN.NU=2^3,
		      MIN.PHI=2^3,
		      THR.NU.PHI=TRUE,
		      thresholdCopynumber=TRUE){
	list(
	     MIN.OBS=MIN.OBS,
	     DF.PRIOR=DF.PRIOR,
	     GT.CONF.THR=GT.CONF.THR,
	     bias.adj=bias.adj,
	     prior.prob=prior.prob,
	     seed=seed,
	     verbose=verbose,
	     PHI.THR=PHI.THR,
	     nHOM.THR=nHOM.THR,
	     MIN.NU=MIN.NU,
	     MIN.PHI=MIN.PHI,
	     THR.NU.PHI=THR.NU.PHI,
	     thresholdCopynumber=thresholdCopynumber)
}

##linear model parameters
lm.parameters <- function(object, cnOptions){
	fD <- fData(object)
	batch <- object$batch
	uplate <- unique(batch)
	parameterNames <- c(paste("tau2A", uplate, sep="_"),
			    paste("tau2B", uplate, sep="_"),
			    paste("sig2A", uplate, sep="_"),
			    paste("sig2B", uplate, sep="_"),
			    paste("nuA", uplate, sep="_"),
			    paste("nuA.se", uplate, sep="_"),			    
			    paste("nuB", uplate, sep="_"),
			    paste("nuB.se", uplate, sep="_"),			    			    
			    paste("phiA", uplate, sep="_"),
			    paste("phiA.se", uplate, sep="_"),			    
			    paste("phiB", uplate, sep="_"),
			    paste("phiB.se", uplate, sep="_"),			    
			    paste("phiAX", uplate, sep="_"),
			    paste("phiBX", uplate, sep="_"),			    
			    paste("corr", uplate, sep="_"),
			    paste("corrA.BB", uplate, sep="_"),
			    paste("corrB.AA", uplate, sep="_"))
	pMatrix <- data.frame(matrix(numeric(0),
				     nrow(object),
				     length(parameterNames)),
				     row.names=featureNames(object))
	colnames(pMatrix) <- parameterNames
	fD2 <- cbind(fD, pMatrix)
	new("AnnotatedDataFrame", data=fD2,
	    varMetadata=data.frame(labelDescription=colnames(fD2),
	    row.names=colnames(fD2)))
}

nonpolymorphic <- function(object, cnOptions, tmp.objects){
	batch <- unique(object$batch)
	CHR <- unique(chromosome(object))
	goodSnps <- function(object, PHI.THR, tmp.objects, nAA.THR, nBB.THR){
		Ns <- tmp.objects[["Ns"]]
		##Ns <- get("Ns", envir)
		flags <- getFlags(object, PHI.THR)
		fewAA <- Ns[, "AA"] < nAA.THR
		fewBB <- Ns[, "BB"] < nBB.THR
		flagsA <- flags | fewAA
		flagsB <- flags | fewBB
		flags <- list(A=flagsA, B=flagsB)
		return(flags)
	}
	nAA.THR <- cnOptions$nHOM.THR
	nBB.THR <- cnOptions$nHOM.THR
	PHI.THR <- cnOptions$PHI.THR
	snpflags <- goodSnps(object, PHI.THR, tmp.objects, nAA.THR, nBB.THR)
	flagsA <- snpflags$A
	flagsB <- snpflags$B
##	if(all(flagsA) | all(flagsB)) stop("all snps are flagged")
	nuA <- getParam(object, "nuA", batch)
	nuB <- getParam(object, "nuB", batch)
	phiA <- getParam(object, "phiA", batch)
	phiB <- getParam(object, "phiB", batch)	
	sns <- sampleNames(object)
	muA <- tmp.objects[["muA"]]
	muB <- tmp.objects[["muB"]]
	A <- A(object)
	B <- B(object)
	CA <- CA(object)
	CB <- CB(object)
	if(CHR == 23){
		phiAX <- getParam(object, "phiAX", batch)
		phiBX <- getParam(object, "phiBX", batch)
	}
	##---------------------------------------------------------------------------
	## Train on unflagged SNPs
	##---------------------------------------------------------------------------
	##Might be best to train using the X chromosome, since for the
	##X phi and nu have been adjusted for cross-hybridization
	##plateInd <- plate == uplate[p]
	##muA <- muA[!flagsA, p, c("A", "AA")]
	##muB <- muB[!flagsB, p, c("B", "BB")]
	muA <- muA[!flagsA, "AA"]
	muB <- muB[!flagsB, "BB"]
	X <- cbind(1, log2(c(muA, muB)))
	Y <- log2(c(phiA[!flagsA], phiB[!flagsB]))
	if(nrow(X) > 5000){
		ix <- sample(1:nrow(X), 5000)
	} else {
		ix <- 1:nrow(X)
	}
	betahat <- solve(crossprod(X[ix, ]), crossprod(X[ix, ], Y[ix]))
	normal <- tmp.objects[["normal"]][!isSnp(object), ]	
	if(CHR == 23){
		##normalNP <- envir[["normalNP"]]
		##normalNP <- normalNP[, plate==uplate[p]]
		##nuT <- envir[["nuT"]]
		##phiT <- envir[["phiT"]]
		
		##cnvs <- envir[["cnvs"]]
                ##loader("cnProbes.rda", pkgname=pkgname, envir=.crlmmPkgEnv)
                ##cnProbes <- get("cnProbes", envir=.crlmmPkgEnv)
		##cnProbes <- cnProbes[match(cnvs, rownames(cnProbes)), ]

		##For build Hg18
		##http://genome.ucsc.edu/cgi-bin/hgGateway
		##pseudo-autosomal regions on X
		##chrX:1-2,709,520 and chrX:154584237-154913754, respectively
		##par:pseudo-autosomal regions
		pseudoAR <- position(object) < 2709520 | (position(object) > 154584237 & position(object) < 154913754)
		##pseudoAR <- cnProbes[, "position"] < 2709520 | (cnProbes[, "position"] > 154584237 & cnProbes[, "position"] < 154913754)
		##in case some of the cnProbes are not annotated
		pseudoAR[is.na(pseudoAR)] <- FALSE
		pseudoAR <- pseudoAR[!isSnp(object)]
		##gender <- envir[["gender"]]
		gender <- object$gender
		obj1 <- object[!isSnp(object), ]
		A.male <- A(obj1[, gender==1])
		mu1 <- rowMedians(A.male, na.rm=TRUE)
		##mu1 <- rowMedians(NP[, gender=="male"], na.rm=TRUE)
		##mu2 <- rowMedians(NP[, gender=="female"], na.rm=TRUE)
		A.female <- A(obj1[, gender==2])
		mu2 <- rowMedians(A.female, na.rm=TRUE)
		mus <- log2(cbind(mu1, mu2))
		X.men <- cbind(1, mus[, 1])
		X.fem <- cbind(1, mus[, 2])
		
		Yhat1 <- as.numeric(X.men %*% betahat)
		Yhat2 <- as.numeric(X.fem %*% betahat)
		phi1 <- 2^(Yhat1)
		phi2 <- 2^(Yhat2)
		nu1 <- 2^(mus[, 1]) - phi1
		nu2 <- 2^(mus[, 2]) - 2*phi2

		if(any(pseudoAR)){
			nu1[pseudoAR] <- 2^(mus[pseudoAR, 1]) - 2*phi1[pseudoAR]
		}
		CT1 <- 1/phi1*(A.male-nu1)
		CT2 <- 1/phi2*(A.female-nu2)
		##CT2 <- 1/phi2*(NP[, gender=="female"]-nu2)
		##CT1 <- matrix(as.integer(100*CT1), nrow(CT1), ncol(CT1))
		##CT2 <- matrix(as.integer(100*CT2), nrow(CT2), ncol(CT2))
		##CT <- envir[["CT"]]
		CA <- CA(obj1)
		CA[, gender==1] <- CT1
		CA[, gender==2] <- CT2
		CA(object)[!isSnp(object), ] <- CA
		##CT[, plate==uplate[p] & gender=="male"] <- CT1
		##CT[, plate==uplate[p] & gender=="female"] <- CT2
		##envir[["CT"]] <- CT

		##only using females to compute the variance
		##normalNP[, gender=="male"] <- NA
		normal[, gender==1] <- NA
		sig2A <- getParam(object, "sig2A", batch)
		normal.f <- normal[, object$gender==2]
		sig2A[!isSnp(object)] <- rowMAD(log2(A.female*normal.f), na.rm=TRUE)^2
		sig2A[!isSnp(object) & is.na(sig2A)] <- median(sig2A[!isSnp(object)], na.rm=TRUE)
		##sig2T[, p] <- rowMAD(log2(NP*normalNP), na.rm=TRUE)^2
		object <- pr(object, "sig2A", batch, sig2A)

		nuA[!isSnp(object)] <- nu2
		phiA[!isSnp(object)] <- phi2
		
		THR.NU.PHI <- cnOptions$THR.NU.PHI
		if(THR.NU.PHI){
			verbose <- cnOptions$verbose
			##Assign values to object
			object <- pr(object, "nuA", batch, nuA)
			object <- pr(object, "phiA", batch, phiA)			
			if(verbose) message("Thresholding nu and phi")
			object <- thresholdModelParams(object, cnOptions)
		} else {
			object <- pr(object, "nuA", batch, nuA)		
			object <- pr(object, "phiA", batch, phiA)
		}
	} else {
		A <- A(object)[!isSnp(object), ]
		mus <- rowMedians(A * normal, na.rm=TRUE)
		crosshyb <- max(median(muA) - median(mus), 0)
		X <- cbind(1, log2(mus+crosshyb))
		logPhiT <- X %*% betahat
		phiA[!isSnp(object)] <- 2^(logPhiT)
		nuA[!isSnp(object)] <- mus-2*phiA[!isSnp(object)]

		THR.NU.PHI <- cnOptions$THR.NU.PHI
		if(THR.NU.PHI){
			verbose <- cnOptions$verbose
			##Assign values to object
			object <- pr(object, "nuA", batch, nuA)
			object <- pr(object, "phiA", batch, phiA)			
			if(verbose) message("Thresholding nu and phi")
			object <- thresholdModelParams(object, cnOptions)
			##reassign values (now thresholded at MIN.NU and MIN.PHI
			nuA <- getParam(object, "nuA", batch)
			phiA <- getParam(object, "phiA", batch)
		}
		CA(object)[!isSnp(object), ] <- 1/phiA[!isSnp(object)]*(A - nuA[!isSnp(object)])
		sig2A <- getParam(object, "sig2A", batch)
		sig2A[!isSnp(object)] <- rowMAD(log2(A*normal), na.rm=TRUE)^2
		object <- pr(object, "sig2A", batch, sig2A)
		##added
		object <- pr(object, "nuA", batch, nuA)
		object <- pr(object, "phiA", batch, phiA)
	}
	return(object)
}

##sufficient statistics on the intensity scale
withinGenotypeMoments <- function(object, cnOptions, tmp.objects){
	normal <- tmp.objects[["normal"]]
	## muA, muB: robust estimates of the within-genotype center (intensity scale)
	muA <- tmp.objects[["muA"]]
	muB <- tmp.objects[["muB"]]
	## vA, vB: robust estimates of the within-genotype variance (intensity scale)
	vA <- tmp.objects[["vA"]]
	vB <- tmp.objects[["vB"]]
	Ns <- tmp.objects[["Ns"]]
	G <- snpCallProbability(object) 
	GT.CONF.THR <- cnOptions$GT.CONF.THR
	CHR <- unique(chromosome(object))

	A <- A(object)
	B <- B(object)
##	highConf <- (1-exp(-confs(object)/1000)) > GT.CONF.THR
	xx <- snpCallProbability(object)
	highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
	##highConf <- confs(object) > GT.CONF.THR
	##highConf <- highConf > GT.CONF.THR
	if(CHR == 23){
		gender <- object$gender
##		gender <- envir[["gender"]]
		IX <- matrix(gender, nrow(G), ncol(G), byrow=TRUE)
##		IX <- IX == "female"
		IX <- IX == 2  ##2=female, 1=male
	} else IX <- matrix(TRUE, nrow(G), ncol(G))
	index <- GT.B <- GT.A <- vector("list", 3)
	names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
	##--------------------------------------------------
	##within-genotype sufficient statistics
	##--------------------------------------------------
	##GT.B <- GT.A <- list()
	snpIndicator <- matrix(isSnp(object), nrow(object), ncol(object)) ##RS: added
	for(j in 1:3){
		GT <- G==j & highConf & IX & snpIndicator
		GT <- GT * normal
		Ns[, j+2] <- rowSums(GT, na.rm=TRUE)				
		GT[GT == FALSE] <- NA
		GT.A[[j]] <- GT*A
		GT.B[[j]] <- GT*B
		index[[j]] <- which(Ns[, j+2] > 0 & isSnp(object)) ##RS: added		
		muA[, j+2] <- rowMedians(GT.A[[j]], na.rm=TRUE)
		muB[, j+2] <- rowMedians(GT.B[[j]], na.rm=TRUE)
		vA[, j+2] <- rowMAD(GT.A[[j]], na.rm=TRUE)
		vB[, j+2] <- rowMAD(GT.B[[j]], na.rm=TRUE)

		##Shrink towards the typical variance
		DF <- Ns[, j+2]-1
		DF[DF < 1] <- 1
		v0A <- median(vA[, j+2], na.rm=TRUE)
		v0B <- median(vB[, j+2], na.rm=TRUE)
		if(v0A == 0) v0A <- NA
		if(v0B == 0) v0B <- NA
		DF.PRIOR <- cnOptions$DF.PRIOR
		vA[, j+2] <- (vA[, j+2]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)
		vA[is.na(vA[, j+2]), j+2] <- v0A
		vB[, j+2] <- (vB[, j+2]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)
		vB[is.na(vB[, j+2]), j+2] <- v0B
	}
	if(CHR == 23){
		k <- 1
		for(j in c(1,3)){
			GT <- G==j & highConf & !IX 
			Ns[, k] <- rowSums(GT)
			GT[GT == FALSE] <- NA
			muA[, k] <- rowMedians(GT*A, na.rm=TRUE)
			muB[, k] <- rowMedians(GT*B, na.rm=TRUE)
			vA[, k] <- rowMAD(GT*A, na.rm=TRUE)
			vB[, k] <- rowMAD(GT*B, na.rm=TRUE)
			
			DF <- Ns[, k]-1
			DF[DF < 1] <- 1
			v0A <- median(vA[, k], na.rm=TRUE)
			v0B <- median(vB[, k], na.rm=TRUE)
			vA[, k] <- (vA[, k]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)
			vA[is.na(vA[, k]), k] <- v0A
			vB[, k] <- (vB[, k]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)
			vB[is.na(vB[, k]), k] <- v0B			
			k <- k+1
		}
	}
	tmp.objects[["Ns"]] <- Ns
	tmp.objects[["vA"]] <- vA
	tmp.objects[["vB"]] <- vB
	tmp.objects[["muA"]] <- muA
	tmp.objects[["muB"]] <- muB
	tmp.objects$index <- index
	tmp.objects$GT.A <- GT.A
	tmp.objects$GT.B <- GT.B
	return(tmp.objects)
}


oneBatch <- function(object, cnOptions, tmp.objects){
	muA <- tmp.objects[["muA"]]
	muB <- tmp.objects[["muB"]]
	Ns <- tmp.objects[["Ns"]]
	CHR <- unique(chromosome(object))
	##---------------------------------------------------------------------------
	## Impute sufficient statistics for unobserved genotypes (plate-specific)
	##---------------------------------------------------------------------------
	MIN.OBS <- cnOptions$MIN.OBS
	index.AA <- which(Ns[, "AA"] >= MIN.OBS)
	index.AB <- which(Ns[, "AB"] >= MIN.OBS)
	index.BB <- which(Ns[, "BB"] >= MIN.OBS)
	correct.orderA <- muA[, "AA"] > muA[, "BB"]
	correct.orderB <- muB[, "BB"] > muB[, "AA"]
	##For chr X, this will ignore the males 
	nobs <- rowSums(Ns[, 3:5] >= MIN.OBS, na.rm=TRUE) == 3
	index.complete <- which(correct.orderA & correct.orderB & nobs) ##be selective here
	size <- min(5000, length(index.complete))
	if(size == 5000) index.complete <- sample(index.complete, 5000)
	if(length(index.complete) < 200){
		stop("fewer than 200 snps pass criteria for predicting the sufficient statistics")
	}
	index <- tmp.objects[["index"]]
	index[[1]] <- which(Ns[, "AA"] == 0 & (Ns[, "AB"] >= MIN.OBS & Ns[, "BB"] >= MIN.OBS))
	index[[2]] <- which(Ns[, "AB"] == 0 & (Ns[, "AA"] >= MIN.OBS & Ns[, "BB"] >= MIN.OBS))
	index[[3]] <- which(Ns[, "BB"] == 0 & (Ns[, "AB"] >= MIN.OBS & Ns[, "AA"] >= MIN.OBS))
	mnA <- muA[, 3:5]
	mnB <- muB[, 3:5]
	for(j in 1:3){
		if(length(index[[j]]) == 0) next()
		X <- cbind(1, mnA[index.complete,  -j, drop=FALSE], mnB[index.complete,  -j, drop=FALSE])
		Y <- cbind(mnA[index.complete, j], mnB[index.complete,  j])
		betahat <- solve(crossprod(X), crossprod(X,Y))
		X <- cbind(1, mnA[index[[j]],  -j, drop=FALSE],  mnB[index[[j]],  -j, drop=FALSE])
		mus <- X %*% betahat
		muA[index[[j]], j+2] <- mus[, 1]
		muB[index[[j]], j+2] <- mus[, 2]
	}
	nobsA <- Ns[, "A"] > 10
	nobsB <- Ns[, "B"] > 10
	notMissing <- !(is.na(muA[, "A"]) | is.na(muA[, "B"]) | is.na(muB[, "A"]) | is.na(muB[, "B"]))
	complete <- list()
	complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
	complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here	
	size <- min(5000, length(complete[[1]]))
	if(size == 5000) complete <- lapply(complete, function(x) sample(x, size))
	if(CHR == 23){
		index <- list()
		index[[1]] <- which(Ns[, "A"] == 0)
		index[[2]] <- which(Ns[, "B"] == 0)
		cols <- 2:1
		for(j in 1:2){
			if(length(index[[j]]) == 0) next()
			X <- cbind(1, muA[complete[[j]], cols[j]], muB[complete[[j]], cols[j]])
			Y <- cbind(muA[complete[[j]], j], muB[complete[[j]], j])
			betahat <- solve(crossprod(X), crossprod(X,Y))
			X <- cbind(1, muA[index[[j]], cols[j]],  muB[index[[j]], cols[j]])
			mus <- X %*% betahat
			muA[index[[j]], j] <- mus[, 1]
			muB[index[[j]], j] <- mus[, 2]
		}
	}
	##missing two genotypes
	noAA <- Ns[, "AA"] < MIN.OBS
	noAB <- Ns[, "AB"] < MIN.OBS
	noBB <- Ns[, "BB"] < MIN.OBS
	index[[1]] <- noAA & noAB
	index[[2]] <- noBB & noAB
	index[[3]] <- noAA & noBB
##	snpflags <- envir[["snpflags"]]
	snpflags <- index[[1]] | index[[2]] | index[[3]]
##	snpflags[, p] <- index[[1]] | index[[2]] | index[[3]]

	##---------------------------------------------------------------------------
	## Two genotype clusters not observed -- would sequence help? (didn't seem to)
	## 1. extract index of complete data
	## 2. Regress  mu_missing ~ sequence + mu_observed
	## 3. solve for nu assuming the median is 2
	##---------------------------------------------------------------------------
	cols <- c(3, 1, 2)
	for(j in 1:3){
		if(sum(index[[j]]) == 0) next()
		k <- cols[j]
		X <- cbind(1, mnA[index.complete, k], mnB[index.complete, k])
		Y <- cbind(mnA[index.complete,  -k],
			   mnB[index.complete,  -k])
		betahat <- solve(crossprod(X), crossprod(X,Y))
		X <- cbind(1, mnA[index[[j]],  k], mnB[index[[j]],  k])
		mus <- X %*% betahat
		muA[index[[j]], -c(1, 2, k+2)] <- mus[, 1:2]
		muB[index[[j]], -c(1, 2, k+2)] <- mus[, 3:4]
	}
	negA <- rowSums(muA < 0) > 0
	negB <- rowSums(muB < 0) > 0	
	snpflags <- snpflags| negA | negB | rowSums(is.na(muA[, 3:5]), na.rm=TRUE) > 0
	tmp.objects$snpflags <- snpflags
	tmp.objects[["muA"]] <- muA
	tmp.objects[["muB"]] <- muB
	tmp.objects[["snpflags"]] <- snpflags
	return(tmp.objects)
}

##Estimate tau2, sigma2, and correlation (updates the object)
locationAndScale <- function(object, cnOptions, tmp.objects){
	DF.PRIOR <- cnOptions$DF.PRIOR
	Ns <- tmp.objects[["Ns"]]
	index <- tmp.objects[["index"]]
	index.AA <- index[[1]]
	index.AB <- index[[2]]
	index.BB <- index[[3]]
	rm(index); gc()

	GT.A <- tmp.objects[["GT.A"]]
	GT.B <- tmp.objects[["GT.B"]]
	AA.A <- GT.A[[1]]
	AB.A <- GT.A[[2]]
	BB.A <- GT.A[[3]]
	
	AA.B <- GT.B[[1]]
	AB.B <- GT.B[[2]]
	BB.B <- GT.B[[3]]
	x <- BB.A[index.BB, ]
	batch <- unique(object$batch)
	tau2A <- getParam(object, "tau2A", batch)
	tau2A[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2	
	DF <- Ns[, "BB"]-1
	DF[DF < 1] <- 1
	med <- median(tau2A, na.rm=TRUE)
	tau2A <- (tau2A * DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
	tau2A[is.na(tau2A) & isSnp(object)] <- med
	object <- pr(object, "tau2A", batch, tau2A)

	sig2B <- getParam(object, "sig2B", batch)	
	x <- BB.B[index.BB, ]
	sig2B[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2	
	med <- median(sig2B, na.rm=TRUE)
	sig2B <- (sig2B * DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
	sig2B[is.na(sig2B) & isSnp(object)] <- med
	object <- pr(object, "sig2B", batch, sig2B)	

	tau2B <- getParam(object, "tau2B", batch)	
	x <- AA.B[index.AA, ]
	tau2B[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2		
	DF <- Ns[, "AA"]-1
	DF[DF < 1] <- 1
	med <- median(tau2B, na.rm=TRUE)
	tau2B <- (tau2B * DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
	tau2B[is.na(tau2B) & isSnp(object)] <- med
	object <- pr(object, "tau2B", batch, tau2B)

	sig2A <- getParam(object, "sig2A", batch)	
	x <- AA.A[index.AA, ]
	sig2A[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2##var(log(IA)|AA)	
	med <- median(sig2A, na.rm=TRUE)
	sig2A <- (sig2A * DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
	sig2A[is.na(sig2A) & isSnp(object)] <- med
	object <- pr(object, "sig2A", batch, sig2A)

	if(length(index.AB) > 0){ ##all homozygous is possible
		corr <- getParam(object, "corr", batch)
		x <- AB.A[index.AB, ]
		y <- AB.B[index.AB, ]
		corr[index.AB] <- rowCors(x, y, na.rm=TRUE)
		corr[corr < 0] <- 0
		DF <- Ns[, "AB"]-1
		DF[DF<1] <- 1
		med <- median(corr, na.rm=TRUE)
		corr <- (corr*DF  +  med * DF.PRIOR)/(DF.PRIOR + DF)
		corr[is.na(corr) & isSnp(object)] <- med
		object <- pr(object, "corr", batch, corr)		
	}
	corrB.AA <- getParam(object, "corrB.AA", batch)
	backgroundB <- AA.B[index.AA, ]
	signalA <- AA.A[index.AA, ]
	corrB.AA[index.AA] <- rowCors(backgroundB, signalA, na.rm=TRUE)
	DF <- Ns[, "AA"]-1
	DF[DF < 1] <- 1
	med <- median(corrB.AA, na.rm=TRUE)
	corrB.AA <- (corrB.AA*DF + med*DF.PRIOR)/(DF.PRIOR + DF)
	corrB.AA[is.na(corrB.AA) & isSnp(object)] <- med
	object <- pr(object, "corrB.AA", batch, corrB.AA)

	corrA.BB <- getParam(object, "corrA.BB", batch)	
	backgroundA <- BB.A[index.BB, ]
	signalB <- BB.B[index.BB, ]
	corrA.BB[index.BB] <- rowCors(backgroundA, signalB, na.rm=TRUE)
	DF <- Ns[, "BB"]-1
	DF[DF < 1] <- 1
	med <- median(corrA.BB, na.rm=TRUE)
	corrA.BB <- (corrA.BB*DF + med*DF.PRIOR)/(DF.PRIOR + DF)
	corrA.BB[is.na(corrA.BB) & isSnp(object)] <- med
	object <- pr(object, "corrA.BB", batch, corrA.BB)
	return(object)
}

coefs <- function(object, cnOptions, tmp.objects){
	batch <- unique(object$batch)
	CHR <- unique(chromosome(object))
	muA <- tmp.objects[["muA"]]
	muB <- tmp.objects[["muB"]]
	vA <- tmp.objects[["vA"]]
	vB <- tmp.objects[["vB"]]
	Ns <- tmp.objects[["Ns"]]
	if(CHR != 23){
		IA <- muA[, 3:5]
		IB <- muB[, 3:5]
		vA <- vA[, 3:5]
		vB <- vB[, 3:5]
		Np <- Ns[, 3:5]
	} else {
		NOHET <- is.na(median(vA[, "AB"], na.rm=TRUE))
		if(NOHET){
			IA <- muA[, -4]
			IB <- muB[, -4]
			vA <- vA[, -4]
			vB <- vB[, -4]
			Np <- Ns[, -4] 
		} else{
			IA <- muA
			IB <- muB
			vA <- vA
			vB <- vB
			Np <- Ns
		}
		
	}
	Np[Np < 1] <- 1
	vA2 <- vA^2/Np
	vB2 <- vB^2/Np
	wA <- sqrt(1/vA2)
	wB <- sqrt(1/vB2)
	YA <- IA*wA
	YB <- IB*wB
	##update lm.coefficients stored in object
	object <- nuphiAllele(object, allele="A", Ystar=YA, W=wA, tmp.objects=tmp.objects, cnOptions=cnOptions)
	object <- nuphiAllele(object, allele="B", Ystar=YB, W=wB, tmp.objects=tmp.objects, cnOptions=cnOptions)	
	##---------------------------------------------------------------------------
	##Estimate crosshyb using X chromosome and sequence information
	##---------------------------------------------------------------------------
	##browser()
	####data(sequences, package="genomewidesnp6Crlmm")
	##snpflags <- envir[["snpflags"]]
	##muA <- envir[["muA"]][, p, 3:5]
	##muB <- envir[["muB"]][, p, 3:5]
	##Y <- envir[["phiAx"]]
	##load("sequences.rda")
	##seqA <- sequences[, "A", ][, 1]
	##seqA <- seqA[match(snps, names(seqA))]
	##X <- cbind(1, sequenceDesignMatrix(seqA))
	##X <- cbind(X, nuA[, p], phiA[, p], nuB[, p], phiB[, p])
	##missing <- rowSums(is.na(X)) > 0
	##betahat <- solve(crossprod(X[!missing, ]), crossprod(X[!missing, ], Y[!missing]))
	return(object)
}

polymorphic <- function(object, cnOptions, tmp.objects){
	batch <- unique(object$batch)
	CHR <- unique(chromosome(object))
	vA <- tmp.objects[["vA"]]
	vB <- tmp.objects[["vB"]]
	Ns <- tmp.objects[["Ns"]]
	
	nuA <- getParam(object, "nuA", batch)
	nuB <- getParam(object, "nuB", batch)
	nuA.se <- getParam(object, "nuA.se", batch)
	nuB.se <- getParam(object, "nuB.se", batch)

	phiA <- getParam(object, "phiA", batch)
	phiB <- getParam(object, "phiB", batch)
	phiA.se <- getParam(object, "phiA.se", batch)
	phiB.se <- getParam(object, "phiB.se", batch)
	A <- A(object)
	B <- B(object)

	NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05
	##---------------------------------------------------------------------------
	## Estimate CA, CB
	##---------------------------------------------------------------------------
	if(CHR == 23){
		phiAX <- getParam(object, "phiAX", batch)  ##nonspecific hybridization coef
		phiBX <- getParam(object, "phiBX", batch)  ##nonspecific hybridization coef			
		phistar <- phiBX/phiA  
		tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB
		copyB <- tmp/(1-phistar*phiAX/phiB)
		copyA <- (A-nuA-phiAX*copyB)/phiA
		CB(object) <- copyB  ## multiplies by 100 and converts to integer 
		CA(object) <- copyA
	} else{
		CA(object) <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))
		CB(object) <- matrix((1/phiB*(B-nuB)), nrow(B), ncol(B))

	}
	return(object)
}

posteriorProbability.snps <- function(object, cnOptions, tmp.objects=list()){
	I <- isSnp(object)
	gender <- object$gender
	CHR <- unique(chromosome(object))
	batch <- unique(object$batch)
	if(CHR == 23){
		phiAX <- getParam(object, "phiAX", batch)[I]
		phiBX <- getParam(object, "phiBX", batch)[I]
	}
	A <- A(object)[I, ]
	B <- B(object)[I, ]
	sig2A <- getParam(object, "sig2A", batch)[I]
	sig2B <- getParam(object, "sig2B", batch)[I]
	tau2A <- getParam(object, "tau2A", batch)[I]
	tau2B <- getParam(object, "tau2B", batch)[I]	
	corrA.BB <- getParam(object, "corrA.BB", batch)[I]
	corrB.AA <- getParam(object, "corrB.AA", batch)[I]
	corr <- getParam(object, "corr", batch)[I]		
	nuA <- getParam(object, "nuA", batch)[I]
	nuB <- getParam(object, "nuB", batch)[I]	
	phiA <- getParam(object, "phiA", batch)[I]
	phiB <- getParam(object, "phiB", batch)[I]
	normal <- tmp.objects[["normal"]][I, ]
	prior.prob <- cnOptions$prior.prob
	emit <- array(NA, dim=c(nrow(A), ncol(A), 10))##SNPs x sample x 'truth'	
	lA <- log2(A)
	lB <- log2(B)	
	X <- cbind(lA, lB)
	counter <- 1##state counter								
	for(CT in 0:3){
		for(CA in 0:CT){
			cat(".")
			CB <- CT-CA
			A.scale <- sqrt(tau2A*(CA==0) + sig2A*(CA > 0))
			B.scale <- sqrt(tau2B*(CA==0) + sig2B*(CA > 0))
			if(CA == 0 & CB == 0) rho <- 0
			if(CA == 0 & CB > 0) rho <- corrA.BB
			if(CA > 0 & CB == 0) rho <- corrB.AA
			if(CA > 0 & CB > 0) rho <- corr
			if(CHR == 23){
				##means <- cbind(suppressWarnings(log2(nuA[, p]+CA*phiA[, p] + CB*phiAx[, p])), suppressWarnings(log2(nuB[, p]+CB*phiB[, p] + CA*phiBx[, p])))
				meanA <- suppressWarnings(log2(nuA+CA*phiA + CB*phiAX))
				meanB <- suppressWarnings(log2(nuB+CB*phiB + CA*phiBX))				
			} else{
				##means <- cbind(suppressWarnings(log2(nuA+CA*phiA)), suppressWarnings(log2(nuB+CB*phiB)))
				meanA <- suppressWarnings(log2(nuA+CA*phiA))
				meanB <- suppressWarnings(log2(nuB+CB*phiB))
				covs <- rho*A.scale*B.scale
				A.scale2 <- A.scale^2
				B.scale2 <- B.scale^2
			}
			Q.x.y <- 1/(1-rho^2)*(((lA - meanA)/A.scale)^2 + ((lB - meanB)/B.scale)^2 - 2*rho*((lA - meanA)*(lB - meanB))/(A.scale*B.scale))
			f.x.y <- 1/(2*pi*A.scale*B.scale*sqrt(1-rho^2))*exp(-0.5*Q.x.y)
			emit[, , counter] <- f.x.y
			counter <- counter+1
		}
	}
	priorProb <- cnOptions$prior.prob
	homDel <- priorProb[1]*emit[, , 1]
	hemDel <- priorProb[2]*emit[, , c(2, 3)] # + priorProb[3]*emit[, c(4, 5, 6)] + priorProb[4]*emit[, c(7:10)]
	norm <- priorProb[3]*emit[, , 4:6]
	amp <- priorProb[4]*emit[, , 7:10]
	##sum over the different combinations within each copy number state
	hemDel <- hemDel[, , 1] + hemDel[, , 2]
	norm <- norm[, , 1] + norm[, , 2] + norm[, , 3]
	amp <- amp[, , 1] + amp[, , 2] + amp[ , , 3] + amp[, , 4]
	total <- homDel + hemDel + norm + amp
	homDel <- homDel/total
	hemDel <- hemDel/total
	norm <- norm/total
	amp <- amp/total
	tmp.objects$posteriorProb <- list(hemDel=hemDel, norm=norm, amp=amp)
	##envir[["posteriorProb"]] <- list(hemDel=hemDel, norm=norm, amp=amp)
	posteriorProb <- array(NA, dim=c(nrow(A), ncol(A), 4))
	posteriorProb[, , 1] <- homDel
	posteriorProb[, , 2] <- hemDel
	posteriorProb[, , 3] <- norm
	posteriorProb[, , 4] <- amp
	return(list(tmp.objects, posteriorProb))
}

biasAdj <- function(object, cnOptions, tmp.objects){
	gender <- object$gender
	CHR <- unique(chromosome(object))
	I <- isSnp(object)
	A <- A(object)[I, ]
	normal <- tmp.objects[["normal"]][I, ]
	results <- posteriorProbability.snps(object=object, cnOptions=cnOptions, tmp.objects=tmp.objects)
	posteriorProb <- results[[2]]
	tmp.objects <- results[[1]]
	mostLikelyState <- apply(posteriorProb, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
	if(CHR == 23){
		##so state index 3 is the most likely state for men and women
		mostLikelyState[, gender==1] <- mostLikelyState[, gender==1] + 1
	}
	proportionSamplesAltered <- rowMeans(mostLikelyState != 3)
	ii <- proportionSamplesAltered < 0.8 & proportionSamplesAltered > 0.01
	##  only exclude observations from one tail, depending on
	##  whether more are up or down
	moreup <- rowSums(mostLikelyState > 3) > rowSums(mostLikelyState < 3) ##3 is normal
	## if equal, which points have a high posterior probability of altered copy number.  drop those.
	NORM <- matrix(FALSE, nrow(A), ncol(A))
	##NORM[proportionSamplesAltered > 0.8, ] <- FALSE
	ratioUp <- posteriorProb[, , 4]/posteriorProb[, , 3]
	NORM[ii & moreup, ] <- ratioUp[moreup & ii] < 1  ##normal more likely
	ratioDown <- posteriorProb[, , 2]/posteriorProb[, , 3]
	NORM[ii & !moreup, ] <- ratioDown[!moreup & ii] < 1  ##normal more likely
	normal <- NORM*normal
	tmp <- tmp.objects[["normal"]]
	tmp[I, ] <- normal
	tmp.objects[["normal"]] <- tmp
	return(tmp.objects)
}


##biasAdjNP <- function(plateIndex, envir, priorProb){
biasAdjNP <- function(object, cnOptions, tmp.objects){
	batch <- unique(object$batch)
	normalNP <- tmp.objects[["normal"]][!isSnp(object), ]
	CHR <- unique(chromosome(object))
	A <- A(object)[!isSnp(object), ]
	sig2A <- getParam(object, "sig2A", batch)
	gender <- object$gender
	##Assume that on the log-scale, that the background variance is the same...
	tau2A <- sig2A
	nuA <- getParam(object, "nuA", batch)
	phiA <- getParam(object, "phiA", batch)
	prior.prob <- cnOptions$prior.prob
	emit <- array(NA, dim=c(nrow(A), ncol(A), 4))##SNPs x sample x 'truth'	
	lT <- log2(A)
	I <- isSnp(object)
	counter <- 1 ##state counter
##	for(CT in 0:3){
##		sds <- sqrt(tau2A[I]*(CT==0) + sig2A[I]*(CT > 0))
##		means <- suppressWarnings(log2(nuA[I]+CT*phiA[I]))
##		tmp <- dnorm(lT, mean=means, sd=sds)
##		emit[, , counter] <- tmp
##		counter <- counter+1
##	}
##	mostLikelyState <- apply(emit, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
	counter <- 1
	for(CT in c(0,1,2,2.5)){
		sds <- sqrt(tau2A[I]*(CT==0) + sig2A[I]*(CT > 0))
		means <- suppressWarnings(log2(nuA[I]+CT*phiA[I]))
		tmp <- dnorm(lT, mean=means, sd=sds)
		emit[, , counter] <- tmp
		counter <- counter+1
	}
	mostLikelyState <- apply(emit, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
	
	if(CHR == 23){
		## the state index for male on chromosome 23  is 2
		## add 1 so that the state index is 3 for 'normal' state
		mostLikelyState[, gender=="male"] <- mostLikelyState[, gender==1] + 1
	}
	tmp3 <- mostLikelyState != 3
	##Those near 1 have NaNs for nu and phi.  this occurs by NaNs in the muA[,, "A"] or muA[, , "B"] for X chromosome
	proportionSamplesAltered <- rowMeans(tmp3)##prop normal
	ii <- proportionSamplesAltered < 0.75
	moreup <- rowSums(mostLikelyState > 3) > rowSums(mostLikelyState < 3)
	notUp <-  mostLikelyState[ii & moreup, ] <= 3
	notDown <- mostLikelyState[ii & !moreup, ] >= 3
	NORM <- matrix(TRUE, nrow(A), ncol(A))
	NORM[ii & moreup, ] <- notUp
	NORM[ii & !moreup, ] <- notDown
	normalNP <- normalNP*NORM

	##flagAltered <- which(proportionSamplesAltered > 0.5)
	##envir[["flagAlteredNP"]] <- flagAltered
	normal <- tmp.objects[["normal"]]
	normal[!isSnp(object), ] <- normalNP
	tmp.objects[["normal"]] <- normal
	return(tmp.objects)
}


getParams <- function(object, batch){
	batch <- unique(object$batch)
	if(length(batch) > 1) stop("batch variable not unique")		
	nuA <- as.numeric(fData(object)[, match(paste("nuA", batch, sep="_"), fvarLabels(object))])
	nuB <- as.numeric(fData(object)[, match(paste("nuB", batch, sep="_"), fvarLabels(object))])	
	phiA <- as.numeric(fData(object)[, match(paste("phiA", batch, sep="_"), fvarLabels(object))])
	phiB <- as.numeric(fData(object)[, match(paste("phiB", batch, sep="_"), fvarLabels(object))])	
	tau2A <- as.numeric(fData(object)[, match(paste("tau2A", batch, sep="_"), fvarLabels(object))])
	tau2B <- as.numeric(fData(object)[, match(paste("tau2B", batch, sep="_"), fvarLabels(object))])
	sig2A <- as.numeric(fData(object)[, match(paste("sig2A", batch, sep="_"), fvarLabels(object))])
	sig2B <- as.numeric(fData(object)[, match(paste("sig2B", batch, sep="_"), fvarLabels(object))])
	corrA.BB <- as.numeric(fData(object)[, match(paste("corrA.BB", batch, sep="_"), fvarLabels(object))])
	corrB.AA <- as.numeric(fData(object)[, match(paste("corrB.AA", batch, sep="_"), fvarLabels(object))])
	corr <- as.numeric(fData(object)[, match(paste("corr", batch, sep="_"), fvarLabels(object))])
	params <- list(nuA=nuA,
		       nuB=nuB,
		       phiA=phiA,
		       phiB=phiB,
		       tau2A=tau2A,
		       tau2B=tau2B,
		       sig2A=sig2A,
		       sig2B=sig2B,
		       corrA.BB=corrA.BB,
		       corrB.AA=corrB.AA,
		       corr=corr)
	return(params)
}


## Constrain nu and phi to positive values
thresholdModelParams <- function(object, cnOptions){
	MIN.NU <- cnOptions$MIN.NU
	MIN.PHI <- cnOptions$MIN.PHI
	batch <- unique(object$batch)
	nuA <- getParam(object, "nuA", batch)
	nuA[nuA < MIN.NU] <- MIN.NU
	object <- pr(object, "nuA", batch, nuA)
	nuB <- getParam(object, "nuB", batch)
	if(!all(is.na(nuB))){
		nuB[nuB < MIN.NU] <- MIN.NU
		object <- pr(object, "nuB", batch, nuB)
	}
	phiA <- getParam(object, "phiA", batch)
	phiA[phiA < MIN.PHI] <- MIN.PHI
	object <- pr(object, "phiA", batch, phiA)
	phiB <- getParam(object, "phiB", batch)
	if(!all(is.na(phiB))){
		phiB[phiB < MIN.PHI] <- MIN.PHI
		object <- pr(object, "phiB", batch, phiB)		
	}
	phiAX <- as.numeric(getParam(object, "phiAX", batch))	
	if(!all(is.na(phiAX))){
		phiAX[phiAX < MIN.PHI] <- MIN.PHI
		object <- pr(object, "phiAX", batch, phiAX)		
	}
	phiBX <- as.numeric(getParam(object, "phiBX", batch))
	if(!all(is.na(phiBX))){
		phiBX[phiBX < MIN.PHI] <- MIN.PHI
		object <- pr(object, "phiBX", batch, phiBX)			
	}
	return(object)
}

##computeCopynumber.SnpSuperSet <- function(object, cnOptions){
####	use.ff <- cnOptions[["use.ff"]]
####	if(!use.ff){
####		object <- as(object, "CrlmmSet")
####	} else	object <- as(object, "CrlmmSetFF")
##	bias.adj <- cnOptions[["bias.adj"]]
##	##must be FALSE to initialize parameters
##	cnOptions[["bias.adj"]] <- FALSE
##	## Add linear model parameters to the CrlmmSet object
##	featureData(object) <- lm.parameters(object, cnOptions)
##	if(!isValidCdfName(annotation(object))) stop(annotation(object), " not supported.")
##	object <- as(object, "CNSet")
##	object <- computeCopynumber.CNSet(object, cnOptions)
##	if(bias.adj==TRUE){## run a second time
##		object <- computeCopynumber.CNSet(object, cnOptions)
##	}
##	return(object)
##}


computeCopynumber.CNSet <- function(object, cnOptions){
	PLATE <- unique(object$batch)
	verbose <- cnOptions$verbose
	tmp.objects <- instantiateObjects(object, cnOptions)
	bias.adj <- cnOptions$bias.adj
	if(bias.adj & ncol(object) <= 15){
		warning(paste("bias.adj is TRUE, but too few samples to perform this step"))
		cnOptions$bias.adj <- bias.adj <- FALSE
	}
	if(bias.adj){
		if(verbose) message("Dropping samples with low posterior prob. of normal copy number (samples dropped is locus-specific)")
		tmp.objects <- biasAdjNP(object, cnOptions, tmp.objects)
		tmp.objects <- biasAdj(object, cnOptions, tmp.objects)
		if(verbose) message("Recomputing location and scale parameters")
	}
	##update tmp.objects
	tmp.objects <- withinGenotypeMoments(object,
					     cnOptions=cnOptions,
					     tmp.objects=tmp.objects)
	object <- locationAndScale(object, cnOptions, tmp.objects)
	tmp.objects <- oneBatch(object,
				cnOptions=cnOptions,
				tmp.objects=tmp.objects)
	##coefs calls nuphiAllele.
	object <- coefs(object, cnOptions, tmp.objects)
	##nuA=getParam(object, "nuA", PLATE)
	THR.NU.PHI <- cnOptions$THR.NU.PHI
	if(THR.NU.PHI){
		verbose <- cnOptions$verbose
		if(verbose) message("Thresholding nu and phi")
		object <- thresholdModelParams(object, cnOptions)
	}		
	if(verbose) message("\nAllele specific copy number")	
	object <- polymorphic(object, cnOptions, tmp.objects)
	if(any(!isSnp(object))){ ## there are nonpolymorphic probes
		if(verbose) message("\nCopy number for nonpolymorphic probes...")	
		object <- nonpolymorphic(object, cnOptions, tmp.objects)
	}
	##---------------------------------------------------------------------------
	##Note: the replacement method multiples by 100
##	CA(object)[, batch==PLATE] <- CA(object)
##	CB(object)[, batch==PLATE] <- CB(object)
	##---------------------------------------------------------------------------
	##update-the plate-specific parameters for copy number
	object <- pr(object, "nuA", PLATE, getParam(object, "nuA", PLATE))
	object <- pr(object, "nuA.se", PLATE, getParam(object, "nuA.se", PLATE))
	object <- pr(object, "nuB", PLATE, getParam(object, "nuB", PLATE))
	object <- pr(object, "nuB.se", PLATE, getParam(object, "nuB.se", PLATE))
	object <- pr(object, "phiA", PLATE, getParam(object, "phiA", PLATE))
	object <- pr(object, "phiA.se", PLATE, getParam(object, "phiA.se", PLATE))
	object <- pr(object, "phiB", PLATE, getParam(object, "phiB", PLATE))
	object <- pr(object, "phiB.se", PLATE, getParam(object, "phiB.se", PLATE))
	object <- pr(object, "tau2A", PLATE, getParam(object, "tau2A", PLATE))
	object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))				
	object <- pr(object, "sig2A", PLATE, getParam(object, "sig2A", PLATE))
	object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))		
	object <- pr(object, "phiAX", PLATE, as.numeric(getParam(object, "phiAX", PLATE)))
	object <- pr(object, "phiBX", PLATE, as.numeric(getParam(object, "phiBX", PLATE)))
	object <- pr(object, "corr", PLATE, getParam(object, "corr", PLATE))
	object <- pr(object, "corrA.BB", PLATE, getParam(object, "corrA.BB", PLATE))
	object <- pr(object, "corrB.AA", PLATE, getParam(object, "corrB.AA", PLATE))
	##object <- object[order(chromosome(object), position(object)), ]
	if(cnOptions[["thresholdCopynumber"]]){
		object <- thresholdCopynumber(object)
	}
	return(object)
}