linearParamElementReplace <- function(obj, elt, value) {
    storage.mode <- storageMode(batchStatistics(obj))
    switch(storage.mode,
           "lockedEnvironment" = {
               aData <- copyEnv(batchStatistics(obj))
               if (is.null(value)) rm(list=elt, envir=aData)
               else aData[[elt]] <- value
               Biobase:::assayDataEnvLock(aData)
               batchStatistics(obj) <- aData
           },
           "environment" = {
               if (is.null(value)) rm(list=elt, envir=batchStatistics(obj))
               else batchStatistics(obj)[[elt]] <- value
           },
           list = batchStatistics(obj)[[elt]] <- value)
    obj
}

## parameters
setMethod("nuA", signature=signature(object="CNSet"), function(object) nu(object, "A"))
setMethod("nuB", signature=signature(object="CNSet"), function(object) nu(object, "B"))
setMethod("phiA", signature=signature(object="CNSet"), function(object) phi(object, "A"))
setMethod("phiB", signature=signature(object="CNSet"), function(object) phi(object, "B"))
setMethod("phiPrimeA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "phiPrimeA")
})
setMethod("phiPrimeB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "phiPrimeB")
})
setMethod("tau2A", signature=signature(object="CNSet"), function(object) tau2(object, "A"))
setMethod("tau2B", signature=signature(object="CNSet"), function(object) tau2(object, "B"))


setMethod("Ns", signature=signature(object="CNSet"),
	   function(object, ...){
		   Ns(batchStatistics(object), ...)
	   })
setMethod("corr", signature=signature(object="CNSet"),
	   function(object, ...){
		   corr(batchStatistics(object), ...)
	   })
setMethod("medians", signature=signature(object="CNSet"),
	   function(object, ...){
		   medians(batchStatistics(object), ...)
	   })
setMethod("mads", signature=signature(object="CNSet"),
	   function(object, ...){
		   mads(batchStatistics(object), ...)
	   })
setMethod("tau2", signature=signature(object="CNSet"),
	   function(object, ...){
		   tau2(batchStatistics(object), ...)
	   })
##---------------------------------------------------------------------------
## Number of samples with Genotype AA, AB, or BB by batch
setMethod("N.AA", signature=signature(object="CNSet"), function(object){
	assayDataElement(batchStatistics(object), "N.AA")
})
setMethod("N.AB", signature=signature(object="CNSet"), function(object){
	assayDataElement(batchStatistics(object), "N.AB")
})
setMethod("N.BB", signature=signature(object="CNSet"), function(object){
	assayDataElement(batchStatistics(object), "N.BB")
})

setReplaceMethod("N.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "N.AA", value)
	  })
setReplaceMethod("N.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "N.AB", value)
	  })
setReplaceMethod("N.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "N.BB", value)
	  })


##---------------------------------------------------------------------------
##  median intensity by genotype cluster for each allele
##should we update the entire matrix rather than one column...
setMethod("medianA.AA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "medianA.AA")
})
setMethod("medianA.AB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "medianA.AB")
})
setMethod("medianA.BB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "medianA.BB")
})
setMethod("medianB.AA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "medianB.AA")
})
setMethod("medianB.AB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "medianB.AB")
})
setMethod("medianB.BB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "medianB.BB")
})
setReplaceMethod("medianA.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "medianA.AA", value)
	  })
setReplaceMethod("medianA.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "medianA.AB", value)
	  })
setReplaceMethod("medianA.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "medianA.BB", value)
	  })
setReplaceMethod("medianB.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "medianB.AA", value)
	  })
setReplaceMethod("medianB.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "medianB.AB", value)
	  })
setReplaceMethod("medianB.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "medianB.BB", value)
	  })

##---------------------------------------------------------------------------
##  mad intensity by genotype cluster for each allele
setMethod("madA.AA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "madA.AA")
})
setMethod("madA.AB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "madA.AB")
})
setMethod("madA.BB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "madA.BB")
})
setMethod("madB.AA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "madB.AA")
})
setMethod("madB.AB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "madB.AB")
})
setMethod("madB.BB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "madB.BB")
})
setReplaceMethod("madA.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "madA.AA", value)
	  })
setReplaceMethod("madA.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "madA.AB", value)
	  })
setReplaceMethod("madA.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "madA.BB", value)
	  })
setReplaceMethod("madB.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "madB.AA", value)
	  })
setReplaceMethod("madB.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "madB.AB", value)
	  })
setReplaceMethod("madB.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "madB.BB", value)
	  })

##---------------------------------------------------------------------------
##  mad of log(intensities) by genotype cluster for each allele
setMethod("tau2A.AA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "tau2A.AA")
})
setMethod("tau2A.BB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "tau2A.BB")
})
setMethod("tau2B.AA", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "tau2B.AA")
})
setMethod("tau2B.BB", signature=signature(object="CNSet"), function(object) {
	assayDataElement(batchStatistics(object), "tau2B.BB")
})
setReplaceMethod("tau2A.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "tau2A.AA", value)
	  })
setReplaceMethod("tau2A.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "tau2A.BB", value)
	  })
setReplaceMethod("tau2B.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "tau2B.AA", value)
	  })
setReplaceMethod("tau2B.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "tau2B.BB", value)
	  })

## correlation of log2A and log2B within each genotype cluster
setMethod("corrAA", signature=signature(object="CNSet"), function(object){
	assayDataElement(batchStatistics(object), "corrAA")
  })
setMethod("corrAB", signature=signature(object="CNSet"), function(object){
	assayDataElement(batchStatistics(object), "corrAB")
  })
setMethod("corrBB", signature=signature(object="CNSet"), function(object){
	assayDataElement(batchStatistics(object), "corrBB")
  })
setReplaceMethod("corrAA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "corrAA", value)
})

setReplaceMethod("corrAB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "corrAB", value)
})

setReplaceMethod("corrBB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "corrBB", value)
})

setReplaceMethod("phiPrimeA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "phiPrimeA", value)
})

setReplaceMethod("phiPrimeB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "phiPrimeB", value)
})

##setMethod("mad.AA", signature=signature(object="CNSet"), function(object) mads(object, "AA"))
##setMethod("mad.AB", signature=signature(object="CNSet"), function(object) mads(object, "AB"))
##setMethod("mad.BB", signature=signature(object="CNSet"), function(object) mads(object, "BB"))
##

##setReplaceMethod("median.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
##	  function(object, value){
##		  linearParamElementReplace(object, "median.AA", value)
##	  })
##setReplaceMethod("median.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
##	  function(object, value){
##		  linearParamElementReplace(object, "median.AB", value)
##	  })
##setReplaceMethod("median.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
##	  function(object, value){
##		  linearParamElementReplace(object, "median.BB", value)
##	  })
##setReplaceMethod("mad.AA", signature=signature(object="CNSet", value="ff_or_matrix"),
##	  function(object, value){
##		  linearParamElementReplace(object, "mad.AA", value)
##	  })
##setReplaceMethod("mad.AB", signature=signature(object="CNSet", value="ff_or_matrix"),
##	  function(object, value){
##		  linearParamElementReplace(object, "mad.AB", value)
##	  })
##setReplaceMethod("mad.BB", signature=signature(object="CNSet", value="ff_or_matrix"),
##	  function(object, value){
##		  linearParamElementReplace(object, "mad.BB", value)
##	  })

setReplaceMethod("nuA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "nuA", value)
	  })

setReplaceMethod("nuB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "nuB", value)
})

setReplaceMethod("phiA", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "phiA", value)
})

setReplaceMethod("phiB", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "phiB", value)
})

setReplaceMethod("tau2A", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "tau2A", value)
})

setReplaceMethod("tau2B", signature=signature(object="CNSet", value="ff_or_matrix"),
	  function(object, value){
		  linearParamElementReplace(object, "tau2B", value)
})



setReplaceMethod("flags", signature=signature(object="CNSet", value="ff_or_matrix"),
		 function(object, value){
			 linearParamElementReplace(object, "flags", value)
})
setReplaceMethod("flags", signature=signature(object="CNSet", value="ff_matrix"),
		 function(object, value){
			 linearParamElementReplace(object, "flags", value)
})


## allele A
##   autosome SNPs
##   autosome NPs
##   chromosome X NPs for women
C1 <- function(object, marker.index, batch.index, sample.index){
##	acn <- vector("list", length(batch.index))
	acn <- matrix(NA, nrow=length(marker.index), ncol=length(sample.index))
	for(k in seq_along(batch.index)){
		l <- batch.index[k]
		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
		bg <- nuA(object)[marker.index, l]
		slope <- phiA(object)[marker.index, l]
		I <- A(object)[marker.index, jj]
		acn[, match(jj, sample.index)] <- 1/slope * (I - bg)
	}
##	if(length(acn) > 1){
##		acn <- do.call("cbind", acn)
##	} else acn <- acn[[1]]
	return(as.matrix(acn))
}

## allele B  (treated allele 'A' for chromosome X NPs)
##   autosome SNPs
##   chromosome X for male nonpolymorphic markers
C2 <- function(object, marker.index, batch.index, sample.index, NP.X=FALSE){
	acn <- matrix(NA, nrow=length(marker.index), ncol=length(sample.index))
	for(k in seq_along(batch.index)){
		l <- batch.index[k]
		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
		bg <- nuB(object)[marker.index, l]
		slope <- phiB(object)[marker.index, l]
		if(!NP.X){
			I <- B(object)[marker.index, jj]
		} else I <- A(object)[marker.index, jj]
		acn[, match(jj, sample.index)] <- 1/slope * (I - bg)
	}
##	if(length(acn) > 1){
##		acn <- do.call("cbind", acn)
##	} else acn <- acn[[1]]
	return(as.matrix(acn))
}

## Chromosome X SNPs
C3 <- function(object, allele, marker.index, batch.index, sample.index){
##	acn <- vector("list", length(batch.index))
	acn <- matrix(NA, nrow=length(marker.index), ncol=length(sample.index))
	for(k in seq_along(batch.index)){
		l <- batch.index[k]
		##j <- which(as.character(batch(object))[sample.index] == batchNames(object)[l])
		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
		phiA2 <- phiPrimeA(object)[marker.index, l]
		phiB2 <- phiPrimeB(object)[marker.index, l]
		phiA <- phiA(object)[marker.index, l]
		phiB <- phiB(object)[marker.index, l]
		nuA <- nuA(object)[marker.index, l]
		nuB <- nuB(object)[marker.index, l]
		IA <- A(object)[marker.index, jj]
		IB <- B(object)[marker.index, jj]
		phistar <- phiB2/phiA
		tmp <- (IB - nuB - phistar*IA + phistar*nuA)/phiB
		CB <- tmp/(1-phistar*phiA2/phiB)
		##CB <- 1/(1-phiA2*phiB2/(phiA*phiB)) * 1/phiB * (IA-nuB-phiB2/phiA*(IA-nuA))
		CA <- (IA-nuA-phiA2*CB)/phiA
		if(allele == "B"){
			acn[, match(jj, sample.index)] <- CB
			##acn[[k]] <- CB
		}
		if(allele == "A"){
			acn[, match(jj, sample.index)] <- (IA-nuA-phiA2*CB)/phiA
		}
		if(allele == "AandB"){
			CA <- tmp/(1-phistar*phiA2/phiB)
			CB <- (IA-nuA-phiA2*CB)/phiA
			acn[, match(jj, sample.index)] <- (IA-nuA-phiA2*CB)/phiA
		}
##		if(allele=="AandB")
##			CA <- tmp/(1-phistar*phiA2/phiB)
##			CB <- (IA-nuA-phiA2*CB)/phiA
##			acn[[k]] <- CA+CB
##		}
	}
##	if(length(acn) > 1){
##		acn <- do.call("cbind", acn)
##	} else acn <- acn[[1]]
	return(as.matrix(acn))
}




ACN <- function(object, allele, i , j){
	if(missing(i) & missing(j)) stop("must specify rows (i) or columns (j)")
	is.ff <- is(calls(object), "ff") | is(calls(object), "ffdf")
	missing.i <- missing(i)
	missing.j <- missing(j)
	if(!missing.i){
		is.ann <- !is.na(chromosome(object)[i])
		is.X <- chromosome(object)[i]==23 & is.ann
		is.auto <- chromosome(object)[i] < 23 & is.ann
		is.snp <- isSnp(object)[i] & is.ann
	} else{
		is.ann <- !is.na(chromosome(object))
		is.X <- chromosome(object)==23 & is.ann
		is.auto <- chromosome(object) < 23 & is.ann
		is.snp <- isSnp(object) & is.ann
		i <- 1:nrow(object)
	}
	## Define batch.index and sample.index
	if(!missing.j) {
		batches <- unique(as.character(batch(object))[j])
		##batches <- as.character(batch(object)[j])
		batch.index <- match(batches, batchNames(object))
	} else {
		batch.index <- seq_along(batchNames(object))
		j <- 1:ncol(object)
	}
	nr <- length(i)
	nc <- length(j)
	acn <- matrix(NA, nr, nc)
	dimnames(acn) <- list(featureNames(object)[i],
			      sampleNames(object)[j])
	if(allele == "A"){
		if(is.ff){
			open(nuA(object))
			open(phiA(object))
			open(A(object))
		}
		## --
		## 4 types of markers for allele A
		##--
		## 1. autosomal SNPs or autosomal NPs
		if(any(is.auto)){
			auto.index <- which(is.auto)
			marker.index <- i[is.auto]
			acn[auto.index, ] <- C1(object, marker.index, batch.index, j)
		}
		if(any(is.X)){
			##2. CHR X SNPs (men and women)
			if(any(is.snp)){
				if(is.ff) {
					open(phiPrimeA(object))
					open(phiPrimeB(object))
					open(phiB(object))
					open(nuB(object))
					open(B(object))
				}
				marker.index <- i[is.X & is.snp]
				acn.index <- which(is.X & is.snp)
				acn[acn.index, ] <- C3(object, allele="A", marker.index, batch.index, j)
				if(is.ff) {
					close(phiPrimeA(object))
					close(phiPrimeB(object))
					close(phiB(object))
					close(nuB(object))
					close(B(object))
				}
			}
			if(any(!is.snp)){
				marker.index <- i[is.X & !is.snp]
				acn.index <- which(is.X & !is.snp)
				acn[acn.index, ] <- NA
				female.index <- j[object$gender[j] == 2]
				## 3. CHR X NPs: women
				if(length(female.index) > 0){
					female.batch.index <- match(unique(as.character(batch(object))[female.index]), batchNames(object))
					jj <- which(object$gender[j] == 2)
					acn[acn.index, jj] <- C1(object, marker.index, female.batch.index, female.index)
				}
##				## 4. CHR X NPs: men
				male.index <- j[object$gender[j] == 1]
				if(length(male.index) > 0){
					if(is.ff){
						open(nuB(object))
						open(phiB(object))
					}
					male.batch.index <- match(unique(as.character(batch(object))[male.index]), batchNames(object))
					jj <- which(object$gender[j] == 1)
					acn[acn.index, jj] <- C2(object, marker.index, male.batch.index, male.index, NP.X=TRUE)
					if(is.ff){
						close(nuB(object))
						close(phiB(object))
					}
				}
			}
		}
		if(is.ff){
			close(nuA(object))
			close(phiA(object))
			close(A(object))
		}
	}
	if(allele == "B"){
		if(is.ff){
			open(nuB(object))
			open(phiB(object))
			open(B(object))
		}
		if(any(!is.snp)){
			acn.index <- which(!is.snp)
			acn[acn.index, ] <- 0
		}
		if(any(is.auto)){
			auto.index <- which(is.auto & is.snp)
			if(length(auto.index) > 0){
				marker.index <- i[auto.index]
				acn[auto.index, ] <- C2(object, marker.index, batch.index, j)
			}
		}
		if(any(is.X)){
			##Chr X SNPs
			if(any(is.snp)){
				if(is.ff) {
					open(phiPrimeA(object))
					open(phiPrimeB(object))
					open(phiA(object))
					open(nuA(object))
					open(A(object))
				}
				marker.index <- i[is.X & is.snp]
				acn.index <- which(is.X & is.snp)
				acn[acn.index, ] <- C3(object, allele="B", marker.index, batch.index, j)
				if(is.ff) {
					close(phiPrimeA(object))
					close(phiPrimeB(object))
					close(phiA(object))
					close(nuA(object))
					close(A(object))
				}
			}
			if(any(!is.snp)){
				acn.index <- which(!is.snp)
				marker.index <- i[!is.snp]
				acn[acn.index, ] <- 0
			}
		}
	}
	if(is.ff){
		close(nuB(object))
		close(phiB(object))
		close(B(object))
	}
	return(acn)
}

setMethod("CA", signature=signature(object="CNSet"),
	  function(object, ...){
		  ca <- ACN(object, allele="A", ...)
		  ca[ca < 0] <- 0
		  ca[ca > 5] <- 5
		  return(ca)
	  })
setMethod("CB", signature=signature(object="CNSet"),
	  function(object, ...) {
		  cb <- ACN(object, allele="B", ...)
		  cb[cb < 0] <- 0
		  cb[cb > 5] <- 5
		  return(cb)
	  })

setMethod("totalCopynumber", signature=signature(object="CNSet"),
	  function(object, ...){
		  ca <- CA(object, ...)
		  cb <- CB(object, ...)
		  return(ca+cb)
	  })

setReplaceMethod("snpCall", c("CNSet", "ff_or_matrix"),
                 function(object, ..., value){
			 assayDataElementReplace(object, "call", value)
		 })
setReplaceMethod("snpCallProbability", c("CNSet", "ff_or_matrix"),
                 function(object, ..., value){
			 assayDataElementReplace(object, "callProbability", value)
		 })