# To use this file:
#  1. Save to a directory on your local hard drive
#  2. Load the file in R using the command:
#     source('c:/path/to/llda.R')
#  3. Run the function, e.g. for the first example (tissue-based):
#     TissueExample<-llda( assay=c(16.5,1,1.65,13), nwell=16, tissue=0.165*10^9, nfp=2, rejl=0.05 )
#     CellExample<-llda( assay=c(0.019,15,0.19,13,1.9,2), nwell=16, cells=19, rejl=0.05 )

### USER REQUIRED VARIABLES ###
# Variable assay: vector of dose levels and number of negative responses
#   i.e. assay<-c( dose1, responses1, dose2, responses2, dose3, responses3, etc)
# Variable nwell: Number of replicate wells
# Variable tissue: Total ng/100ul of tissue [Tissue mode only]
# Variable nfp: Number of footpads [Tissue mode only]
# Variable cells: Total number of cells (x10^6) [Cell mode only]
# Variable rejl: P-value at which to reject hypothesis (default: 0.05)
################################

llda<-function( assay, nwell, tissue=NA, nfp=NA, cells=NA, rejl=0.05 ) {

	# Are we operating in tissue-weight or cell mode?
	UseCells<-NA
	if ( is.na(tissue) && is.na(nfp) && !is.na(cells) ) {
		UseCells<-TRUE
	} else if ( !is.na(tissue) && !is.na(nfp) && is.na(cells) ) {
		UseCells<-FALSE
	} else {
		stop("Unable to determine whether to run in tissue-weight or cells mode.\n - Tissue-weight mode requires values for the tissue and nfp parameters.\n - Cells mode requires values for the cells parameter.")
	}

	# Check input parameter: assay
	if (!is.matrix(assay)) {
		# Turn into a matrix
		if ( length(assay)%%2==0 ) {
			assay<-matrix( assay, byrow=TRUE, ncol=2, nrow=length(assay)/2 )
		} else {
			# Vector length is not divisible by 2
			stop("Assay input is not divisible by two. Correct format is either a 2-column matrix or a vector with alternating doses and negative responses.")
		}
	}
	# Have assay as a matrix now, Remove rows where all wells had all-negative or all-positive responses
	assay<-assay[assay[,2] != nwell,]
	assay<-assay[assay[,2] != 0,]

	if ( UseCells ) {
		# Check input parameter: cells
		if (is.na(cells) || cells<=0) {
			stop("Number of cells (cells) must be greater than 0.")
		}
	} else {
		# Check input parameter: tissue
		if (is.na(tissue) || tissue<=0) {
			stop("Tissue input must be greater than 0, and expressed in terms of ng per 100ul.")
		}

		# Check input parameter: nfp
		if (is.na(nfp) || nfp<=0) {
			stop("Number of footpads (nfp) must be greater than 0.")
		}
	}

	# Check input parameter: nwell
	if (is.na(nwell) || nwell<=0) {
		stop("Number of wells (nwell) must be greater than 0.")
	} else if ( nwell < max(assay[,2]) ) {
		stop("Number of wells (nwell) must be greater or equal to the number of negative responses in the assay.")
	}

	# Check input parameter: rejl
	if (is.na(rejl) || rejl<=0 || rejl>1  ) {
		stop("Rejection level (rejl) must be greater than 0 and less than 1.")
	}
	
	# Block to calculate values needed for computations
	# Variable dd: number of rows in assay
	dd<-nrow(assay)
	# Variable l: List of dose levels (from first column of assay)
	l<-assay[,1]
	# Variable r: List of negative responses (from second column of assay)
	r<-assay[,2]
	# Variable n: creates a new matrix with each 
	# Use of repeat: repeat(M, r, c) matrix that repeats block M, r times in row and c times in column 
	n<-matrix(nwell, nrow=dd, ncol=1 )
	# Variable p: Probability of a negative responses
	p<-r/n
	
	##########################################################################################################################
	## Block for the Assay validity Test
	## Note: Correct according to the example
	wl<-(n * p * (log(p)^2)) / ( 1 - p )
	w2<-(n * p * l^2 ) / ( 1 - p )
	# Variable w: concat by columns w1 and w2
	w<-cbind(wl, w2)
	xl<-l
	x2<-1 / l
	# Variable xx: concat by columns x1 and x2
	xx<-cbind(xl, x2)
	wxx<-w * xx
	
	sumwxx<-apply( wxx, 2, sum )
	sumw<-apply( w, 2, sum )
	xbar<-sumwxx / sumw
	
	x<-xx - matrix(xbar, nrow=dd, ncol=2, byrow=TRUE )
	
	phihat<- -(log(p) / l)
	yy<-cbind( log(phihat), phihat )
	wyy<-w * yy
	sumwyy<-apply( wyy, 2, sum )
	ybar<-sumwyy / sumw

	y<-yy - matrix(ybar, nrow=dd, ncol=2, byrow=TRUE )
	bhat<-apply(w * x * y, 2, sum) / apply(w * x^2, 2, sum)
	vbhat<-abs( ( apply(w * y^2, 2, sum) - apply(w * x * y, 2, sum) ) / ( dd - apply(w * x^2, 2, sum) ) )
	# Variable chi: Chi squared values of AVT
	chi<-bhat * apply(w * x * y, 2, sum)
	
	
	##########################################################################################################################
	## Block for Weighted Least Squares Estimates
	w<-w2
	sumw<-apply( w, 2, sum )
	x<-x2
	y<-y[,2]
	phiwm<-apply(w * phihat, 2, sum) / sumw
	vphiwm<-1/sumw
	phiest<-phiwm
	
	##########################################################################################################################
	## Block for Estimate Validity Test
	## Note: this function uses many global variables
	validest<-function(phiest) {
	
		p<-exp( -1 * phiest * l )
		w1<-( n * p * log(p)^2 ) / (1-p)
		w2<-( n * p * l^2 ) / (1-p)
		w<-cbind( w1, w2 ) 
		x1<-l
		x2<-1/l
		xx<-cbind( x1, x2 )
		wxx<-w * xx

		sumwxx<-apply( wxx, 2, sum )
		sumw<-apply( w, 2, sum )
		xbar<-sumwxx / sumw
	
		x<-xx - matrix(xbar, nrow=dd, ncol=2, byrow=TRUE )
		# yy uses the old phihat, not a new one calculated from the new p.
		yy<-cbind( log(phihat), phihat )
		wyy<-w * yy
		sumwyy<-apply( wyy, 2, sum )
		ybar<-sumwyy / sumw
		y<-yy - matrix(ybar, nrow=dd, ncol=2, byrow=TRUE )
		bhat<-apply(w * x * y, 2, sum) / apply(w * x^2, 2, sum)
	
		# Variable: vbhat (one of the return values)
		vbhat<-abs( ( apply( w * y^2, 2, sum ) - apply(w * x * y, 2, sum) ) / (dd - apply(w * x^2, 2, sum) ) )
		# Variable: chis1 (one of the return values)
		chis1<- bhat * apply(w * x * y, 2, sum)
		chis1<-max(chis1)
	
		return( list(vbhat=vbhat, chis1=chis1) )
	
	}
	
	##########################################################################################################################
	## Block for Chi squared test
	chisq<-function(phimc) {
	
		expon<-phimc * l
		negxpon<-exp(-1 * expon)
	
		return( list(
			negxpon = negxpon,
			negxpon2= exp(-2 * expon),
			negxpon3= exp(-3 * expon),
			posxpon = exp(expon),
			denom    = n * (1 - negxpon)
		) )
	
	}
	
	##########################################################################################################################
	## Block for First Derivative
	first<-function( x ) {
		# x is a list containing negxpon and posxpon
	
		numerf<- n * l * x$negxpon * ( 2 * r - n ) + r^2 * l * ( x$posxpon - 2 )
		jacob<-apply(numerf / (n * (1-x$negxpon)^2), 2, sum)
		return( list(jacob=jacob, numerf=numerf) )
	
	}
	
	##########################################################################################################################
	## Block for Second Derivative
	second<-function(x) {
		# x is a list containing negxpon, negxpon2, negxpon3 and posxpon
		numers<- n * l^2 * (n - 2 * r) * ( x$negxpon - x$negxpon3 ) + r^2 * l^2 * (x$posxpon - 4 + 7 * x$negxpon - 4 * x$negxpon2)
		hess<-apply( numers / (n * (1 - x$negxpon)^4), 2, sum )
	
		return( list(hess=hess, numers=numers) )
	
	}
	
	##########################################################################################################################
	## Block for Newton's Iterative Method
	newton<-function(phiwm) {
	
		# Variable diff: Convergence value
		diff<--0.000001
		phimc<-phiwm
		iter<-0
		check<-FALSE
		chisq_results<-chisq( phimc )
		while ( !check ) {
			iter<-iter+1
			first_results<-first( chisq_results )
			second_results<-second( chisq_results )
			delta<-0.001 * ( first_results$jacob / second_results$hess )
			converge<-phimc - delta
			phimc<-converge
			check<- delta > diff
			chisq_results<-chisq( phimc )
		}
	
		return( list(
			phimc=phimc,
			hess=second_results$hess,
			iter=iter
		))
	
	}
	#rm(list=ls())
	
	##########################################################################################################################
	## Block for main()
	
	## Run Newton's Iterative Method
	newton_results<-newton(phiwm)
	phimc<-newton_results$phimc
	
	## Run Validity Test with phiest from Newton's Iterative method
	validest_results2<-validest(phimc) # Will return vbhat and chis1
	
	## Calculate stats
	# Variable varmc: Estimated variance
	varmc<-2 * 1/newton_results$hess
	# Variable u95 and l95: the 95% upper and lower confidence levels # TODO: Is this from the beginning, or from newton?
	u95<-phimc + 1.96 * sqrt( varmc )
	l95<-phimc - 1.96 * sqrt( varmc )
	# Variable chimc: Chi squared EVT
	chimc<-validest_results2$chis1
	# Variable pv: P-values of AVT
	pv<-1-pchisq( chi, 1 )
	# Variable pvmc: P-value of EVT
	pvmc<-1-pchisq( chimc, 1 )
	# Variable parasite: Number of parasites per lesion
	if ( UseCells ) {
		parasite<-(cells * phimc)
	} else {
		parasite<-tissue/nfp * phimc
	}
	# Variable inv: 1/freq
	inv<-1/phimc
	# Variable iter: iterations to converge
	iter<-newton_results$iter
	
	##########################################################################################################################
	## Is assay valid?
	isAssayValid<-paste("Assay is not valid because the null hypothesis was rejected at the",rejl,"alpha level.")
	if ( all(pv > rejl) ) {
		isAssayValid<-paste("Assay is valid as null hypothesis was not rejected at the",rejl,"alpha level.")
	}
	## Is estimate valid?
	isEstimateValid<-paste("Estimate is not valid because the null hypothesis was rejected at the",rejl,"alpha level.")
	if ( pvmc > rejl ) {
		isEstimateValid<-paste("Estimate is valid as null hypothesis was not rejected at the",rejl,"alpha level.")
	}
	
	##########################################################################################################################
	## Package results
	LModel<-list(
		AVT=list(
			chi=chi,
			pv=pv
		),
		EVT=c(
			estfreq=phimc,
			chi=chimc,
			pv=pvmc,
			estvar=varmc,
			u95=u95,
			l95=l95,
			freqinv=inv,
			parasites=parasite,
			tissue=tissue,
			cells=cells,
			iterations=iter
		)
	)
	
	
	##########################################################################################################################
	## Print results
	printLModel<-function() {
		cat( "\n\t  Minimum Chi-Squared Estimation With Newton's Method\n" )
		cat( "\n-------------------------       Summary       -------------------------\n\n")

		# Summary text depends on whether the assay and estimate was valid.
		if ( all(pv <= rejl) ) {
			# Data does not match our assumptions
			cat( "    Model does not fit the data, and so cannot be used to estimate the frequency or the number of parasites.\n" )
		} else if ( pvmc <= rejl ) {
			cat( "    Model fits the data, but the estimate is not valid, and so cannot be used to estimate the frequency or number of parasites.\n" )
		} else {
			# Assay and Estimate OK, summary depends if we use cells or not
			if ( UseCells ) {
				cat( "    Estimated number of parasites is",signif(parasite,3),"and there are",signif((phimc/10^6),3),"\n    parasites per cell.\n    Models fits the data and the estimated frequency is valid.\n" )
			} else {
				cat( "    Estimated number of parasites is",signif(parasite,3),"and there are",signif(phimc,3),       "\n    parasites per nanogram of tissue.\n    Models fits the data and the estimated frequency is valid.\n" )
			}
		}

		cat( "\n------------------------- Assay Validity Test -------------------------\n\n")
		cat( "\t\t\t\t\tModel1\t\tModel2\n" )
		cat( "    chi-squared AVT\t\t",chi,"\n", sep="\t" )
		cat( "    p-values of AVT\t\t",pv ,"\n", sep="\t" )
		cat( "\n------------------------ Estimate Validity Test -----------------------\n\n")
		if ( UseCells ) {
			cat( "    estimated frequency of parasites per cell",(phimc/10^6),"\n", sep="\t" )
			cat( "    estimated parasites per cell\t",parasite,"\n", sep="\t" )
		} else {
			cat( "    estimated frequency of parasites per lesion",phimc,"\n", sep="\t" )
			cat( "    estimated parasites per lesion\t",parasite,"\n", sep="\t" )
		}
		cat( "    chi-squared EVT\t\t\t",chimc,"\n", sep="\t" )
		cat( "    p-value of EVT\t\t\t",pvmc,"\n", sep="\t" )
		cat( "    estimated variance\t\t\t",varmc,"\n", sep="\t" )
		cat( "    95% upper confidence level\t\t",u95,"\n", sep="\t" )
		cat( "    95% lower confidence level\t\t",l95,"\n", sep="\t" )
		cat( "    1/frequency\t\t\t\t",inv,"\n", sep="\t" )
		cat( "    nanograms of tissue\t\t\t",tissue,"\n", sep="\t" )
		if ( UseCells ) {
			cat( "    total number of cells (10^6)\t",cells,"\n", sep="\t" )
		}
		cat( "    number of iterations to converge\t",iter,"\n", sep="\t" )
		cat( "\n" )
	}
	printLModel()

	return( LModel )
}