"vitality.fit"<-
function(time,sdata,rc.data=F,se=F,gfit=F,datatype="CUM",ttol=.000001,init.params=F,
												 pplot=T,tlab="days",lplot=T,cplot=F,silent=F)
#
#  
#
#     Vitality based survival model: parameter fitting routine:    VERSION: 2/17/2004
#
# REQUIRED PARAMETERS: 
#     time - time component of data: time from experiment start.  Time should 
#	     start after the imposition of a stressor is completed.
#     sdata - survival or mortality data.  The default expects cumulative 
#	     survival fraction.  If providing incremental mortality fraction 
#          instead, use option: datatype="INC".
#          The default also expects the data to represent full mortality.   
#          Otherwise, use option: rc.data=T to indicate right censored data.
#
# OPTIONAL PARAMETERS:
#     rc.data =T  - specifies Right Censored data.   If the data does not 
#          represent full mortality, it is probably right censored.  The default 
#          is rc.data=F.  A third option is rc.data="TF".  Use this case to add
#          a near-term zero survival point to data which displays nearly full 
#          mortality ( <.01 survival at end).  If rc.data=F but the data does
#	     not show full mortality, rc.data="TF" will be 
#          invoked automatically. 
#     se =<population>  calculates the standard errors for the MLE parameters.  
#          Default is se=F.  The initial study population is necessary for 
#          computing these standard errors. 
#     gfit =<population>  provides a Pearson C type test for goodness of fit.  
#	     Default is gfit=F. The initial study population is necessary for 
#	     computing goodness of fit. 
#     datatype ="CUM" -cumulative survival fraction data- is the default.   
#          Other option: datatype="INC" - for incremental mortality fraction 
#          data.  ttol (stopping criteria tolerence.)  Default is .000001 . 
#          specify as ttol=.0001.
#          If one of the liklihood plots (esp. for "k") does not look optimal, 
#	     try decreasing ttol.   If the program crashes, try increasing ttol.
#     init.params =F  has the routine choose initial parameter estimates for 
#          r,s,k  (default: =F). If you wish to specify initial param values 
#          rather than have the routine choose them, specify 
#          init.params=c(r,s,k) in that order (eg. init.params=c(.1,.02,.003)).
#     pplot =T provides plots of cumulative survival and incremental mortality - 
#	     for both data and fitted curves (default: =T).  pplot=F provides no 
#	     plotting.  A third option:  pplot=n  (n>=1) extends the time axis of 
#	     the fitting plots 	(beyond the max time in data).  For example: 
#          pplot=1.2 extends the time axis by 20%.  (Note:  the incremental 
#          mortality plot is a continuous representation of the appropriately-
#          binned histogram of incremental mortalities.)
#     tlab ="<time units>" specifies units for x-axis of plots.  Default is 
#          tlab="days".
#     lplot =T provides likelihood function plotting (default =T).  
#          Note:  these plots are not "likelihood profiles" in that while one 
#          parameter is varied, the others are held fixed, rather than 
#          re-optimized. (must also have pplot=T.)
#     cplot =T provides a likelihood contour plot for a range of r and s values 
#          (can be slow so default is F).   Must also have lplot=T (and pplot=T) 
#          to get contour plots.
#     silent =T stops all print and plot options (still get most warning and all 
#          error messages) Default is F.  A third option, silent="verbose" also 
#          enables the trace setting in the ms (minimum sum) S-Plus routine. 
#
# RETURN:
#     vector of final MLE r,s,k parameter estimates.
#     standard errors of MLE parameter estimates (if se=<population> is
#     specified).
#
{

# --check/prepare flags
	ttrace _F
	if (silent != F) {
		if (silent == T) {
			pplot_ttrace_lplot_F
		}
		else {  # "verbose" mode
			ttrace _T
			silent _F
		}
	}
	if (se == T) {
		se _F
		warning("\n INPUT ERROR: For standard error calculation, you MUST set: se = population size") 
	}
	if (gfit == T) {
		gfit _F
		warning("\n INPUT ERROR: For goodness-of-fit calculation, you MUST set: gfit = population size") 
	}
	
#  --Check/prepare Data---
	dTmp _dataPrep(time,sdata,datatype,rc.data)
	time _dTmp$time
	sfract _dTmp$sfract
	x1 _dTmp$x1
	x2 _dTmp$x2
	Ni _dTmp$Ni
	rc.data _dTmp$rc.data
	
#  --Produce initial parameter values---
	if(length(init.params)==1) {  # use init param estimation scheme
		rsk <-rsk.init(time,sfract)   
	}
	else {                        # use user specified init params
		rsk <-init.params
	}
	if(rsk[1] == -1) {
		stop
	}
	if (silent == F) {
	  print(c("Initial r s k:",rsk))
	}


#  --create dataframe for ms---
	dtfm <- data.frame(x1=x1,x2=x2,Ni=Ni)
	param(dtfm,"r")    
	param(dtfm,"s")   
	param(dtfm,"k")

#  --run MLE fitting routine (ms "minimum sum")---    
   vfit.mle <- ms( ~logLikelihood(r,s,k,x1,x2,Ni),data=dtfm,start=list(r=rsk[1],s=rsk[2],k=rsk[3]),
                   control=list(rel.tolerance=ttol),trace=ttrace)

	if (silent ==F) {print(summary(vfit.mle))}
	
#  --if k < 0, rerun with k=0---
	if (vfit.mle$parameters[3] < 0) {
		k.final _0 
		names(k.final) _"k"
		warning("WARNING: k<0 on initial run.  Trying again with k=0.")
		
		# --create dataframe 
		dtfm <- data.frame(x1=x1,x2=x2,Ni=Ni)   
		param(dtfm,"r")    
		param(dtfm,"s")
		   
		# --run ME
    	vfit.mle <- ms( ~logLikelihood(r,s,0,x1,x2,Ni),data=dtfm,start=list(r=rsk[1],s=rsk[2]),
                   control=list(rel.tolerance=ttol),trace=ttrace)

	   if (silent ==F) {print(summary(vfit.mle))}

		#lplot_F  # k=0 would crash the likelihood plot routine	
	}
	else{
		k.final_ vfit.mle$parameters[3]
	}

# --save final param estimates---
	r.final_ vfit.mle$parameters[1]
	s.final_ vfit.mle$parameters[2]   
	mlv_ vfit.mle$value

#  ==end MLE fitting===


#  --compute standard errors---
	if (se != F) {
		s.e. _stdErr(r.final[1], s.final[1], k.final[1], x1, x2, Ni, se)
	}

	
#  --plotting and goodness of fit---
if (pplot != F || gfit != F) {
	plotting(r.final,s.final,k.final,mlv,time,sfract,x1,x2,Ni,pplot,tlab,lplot,cplot,gfit,rc.data)
}

# ............................................................................................
			
# --return final param values---
	sd_5   #significant digits of output
	if(se != F ) {
		params_c(r.final,s.final,k.final)
		out_signif(cbind( params,s.e.),sd)
		return(out)
	}
	else {
		return(signif(c(r.final,s.final,k.final),sd))
	}
}

# =====================================================================================================

"dataPrep"<-
function(time,sdata,datatype,rc.data)
#  
#  Function to deal with NAs, right truncated data and datatype cumulative survival or
#  incremental motality. 
{
		#check for and remove NAs from data
	if (any(naT_is.na(time))) {
		time   _time[!naT]
		sdata _sdata[!naT]
		warning(message="WARNING:  NAs found in data and removed.")
	}
	if (any(naT_is.na(sdata))) {
		time   _time[!naT]
		sdata _sdata[!naT]
		warning(message="WARNING:  NAs found in data and removed.")
	}
	
	if(length(time) < 3) {stop(message="ERROR:  not enough data.")}
	
	if(!all(0<=sdata) || !all(sdata<=1)) {
		stop("ERROR:  survival fraction data outside range.")
	}
# end data checking	

	maxx2 _max(time)  #for right-censored data and for plotting 
# === check data type (CUMulative or INCremental.  If CUM, create INC ===
if (datatype == "CUM") {
	#====== survival assumed =1 at time 0 ===
	if (time[1] > 0) {
		time _c(0,time)
		sdata _c(1,sdata)
	}
	else {
		if (sdata[1] < 1) {
			sdata_sdata/sdata[1]
			warning(message="Initial survival < 1.  Data scaled so that initial survival =1.")			
		}
	}
	#------------------------
	sfract _sdata
	len_length(time)
	#  ...set up data for MLE fitting of incremental survivorship...	

	# --------- right censored data?
	if (rc.data != T) {
		if (rc.data == F) {
			#check if final sdata indicates full mort
			if (sfract[len] != 0) {
				warning("WARNING: Survival data may be right censored...")
				rc.data_"TF"
			}
			else {
				#standard setup
				x1 _time[1:(len-1)]
				x2 _time[2:len]
				sfract1 _sfract[1:(len-1)]
				sfract2 _sfract[2:len]
			}
		}
		if (rc.data == "TF") {
			#setup: add zero sruv and short time step  ("TF" option)
			x1 _time
			x2 _c(time[2:len],2*time[len]-time[len-1])
			sfract1 _sfract
			sfract2 _c(sfract[2:len],0)
		}
		
	}
	else {                      #if rc.data == T
		x1 _time
		x2 _c(time[2:len],10*maxx2)         #testing...
		sfract1 _sfract
		sfract2 _c(sfract[2:len],0)
	}
	# -----------------  end of dealing with right censored data options	
         
	Ni _sfract1-sfract2     #incrimental survival fraction
	
	#  ...end conversion of cumulative survivorship to incremental mortality
}
else {
	if (datatype == "INC") {
		lent_length(time)
		#should be no t=0 data.  Eliminate if necessary.
		if(time[1] == 0) {
			time  _time[2:lent]
			sdata _sdata[2:lent]
			lent  _length(time)
		}
		#check for right censored data
		if (rc.data != T) {
			if (rc.data == F) {
				if(sum(sdata) < 1) {
					rc.data _"TF"
					warning("WARNING: Survival data may be right censored...")
				}
				else {  
					#standard setup
					Ni _sdata
					x1 _c(0,time)[1:lent]
					x2 _time
				}
			}
			if (rc.data == "TF") {
				#setup: add zero sruv and short time step  ("TF" option)
				Ni _c(sdata,1-sum(sdata))
				x1 _c(0,time)
				x2 _c(time,2*time[lent]-time[lent-1])  #final time interval assumed same as prev.
			}
			
		}
		else {       #if rc.data == T
			Ni _c(sdata,1-sum(sdata))
			x1 _c(0,time)
			x2 _c(time,10*maxx2)    #final time interval large
		}
		
		# Build cumulative data	
	  	#sfract _1.0      # initial survival fraction assumed 1
	    sfract_NULL
	  	for (i in 1:length(Ni)) {
			sfract _c(sfract,1-sum(Ni[1:i]))
	  	}
	time_c(0,time)
	len_length(time)   #check 
	print(c(length(sfract),length(time),length(Ni),length(x1),length(x2),lent))
	
	
		if(!all(-.0001<=sfract<=1)) {    # the -.0001 rather than 0 is to allow some rounding error
		stop("ERROR:  survival fraction outside range.")
		}
   	}
   else {
		stop("ERROR:  bad datatype specification")
	}
}

return(time,sfract,x1,x2,Ni,rc.data)
}

#======================================================================================================

"rsk.init"<-
function(tt,sf)
#
#  Function to find initial r,s,k estimates.    
# tt is time,   
# sf is survival fraction.
{
ii _indexFinder(sf,0.5)
if (ii == -1) {
	warning("ERROR: no survival fraction data below the .5 level.   Can not use the initial r s k estimator.  You must supply initial r s k estimates")
	return(-1)
}
		slope _ (sf[ii]-sf[ii-1]) /(tt[ii]-tt[ii-1])
		t50 _  tt[ii] + (0.5-sf[ii])/slope
		nslope _ slope*t50

rp _r.s.slope[,1]
sp _r.s.slope[,2]
slope.p _r.s.slope[,3]

sze _length(slope.p)
if (slope.p[sze] < nslope && nslope < slope.p[1])  { # check if normalized slope (nslope) is on the chart.
	for (i in 2:sze) {       
		if ( (slope.p[i-1] - nslope)*(slope.p[i] - nslope) < 0.0) {
			rri _( rp[i-1] + (rp[i]-rp[i-1])*(nslope-slope.p[i-1])/(slope.p[i]-slope.p[i-1]) )/t50
			ssi _( sp[i-1] + (sp[i]-sp[i-1])*(nslope-slope.p[i-1])/(slope.p[i]-slope.p[i-1]) )/sqrt(t50)
			break
		}
	}
}
else {    
	if (nslope <= slope.p[sze]) {
		rri _rp[sze]/t50
		ssi _sp[sze]/sqrt(t50)
	}
	else {
		rri _rp[1]/t50
		ssi _sp[1]/sqrt(t50)
	}
}

ssi _ssi/1.1 # ssi was consistently overestimated above.

#  --estimate initial k---

#use rri,ssi and a data point (tt[ii],sf[ii]) to solve for kki 
#      ..using the actual survival function.
ii _indexFinder(sf,.94)-1
if (ii <= 1) {     # In case no data points between sruv=1 and surv=.94.
	ii_ 2
	warning(message="WARNING: Initial time step may be too long.")
}
kki _-(1/tt[ii])*log(sf[ii]/SurvFn(tt[ii],rri,ssi,0))
if (kki <= 0) {
	kki _ (1.0 - sf[ii])/tt[ii] 
}

	
return(c(rri,ssi,kki))
}

#=======================================================================================================

# Script for setting r.s.slope - the data frame used by function rsk.init
# to produce initial r,s estimates from an estimated slope of the survival 
# curve at the inflection point.
c1_c(0.000,0.100,0.200,0.300,0.400,0.500,0.600,0.700,0.800,0.900,0.930,0.950,0.960,0.970,
     0.980,0.990,0.991,0.992,0.993,0.994,0.995,0.996,0.997,0.998)
c2_c(1.48260200,1.40208800,1.31746800,1.22796200,1.13249500,1.02951700,0.91664920,0.78987300,0.64132740,
0.45059940,0.37620130,0.31747940,0.28374720,0.24554230,0.20032630,0.14153800,0.13426380,0.12657460,
0.11839010,0.10959890,0.10004140,0.08947239,0.07747895,0.06325606)
c3_c(-0.2143371,-0.2315594,-0.2518273,-0.2761605,-0.3061416,-0.3443959,-0.3956923,-0.4699249,-0.5925325,
-0.8638228,-1.0422494,-1.2411051,-1.3920767,-1.6126579,-1.9815621,-2.8115970,-2.9646636,-3.1455461,
-3.3638416,-3.6345703,-3.9827944,-4.4543775,-5.1451825,-6.3036320)
r.s.slope_data.frame(c1,c2,c3)
dimnames(r.s.slope)[[2]]_c("r","s","slope")
rm(c1)
rm(c2)
rm(c3)

#=========================================================================================================

"SurvFn"<-
function(xx,r,s,k) 
 #  The cumulative survival distribution function.
{
	# pnorm is: cumulative prob for the Normal Dist.
	tmp1 <- (sqrt(1/xx) * (1 - xx * r))/s    #  xx=0 is ok.  pnorm(+-Inf) is defined
	tmp2 <- (sqrt(1/xx) * (1 + xx * r))/s

	# --safeguard if exponent gets too large.---
	tmp3 <- 2*r/(s*s)
	
	if (tmp3 >250) {   
		q _tmp3/250 
		
		if (tmp3 >1500) {
			q _tmp3/500
		}
		
		valueFF <-(1.-(pnorm(-tmp1) + (exp(tmp3/q) *pnorm(-tmp2)^(1/q))^(q)))*exp(-k*xx) 
	}		    
	else {
		valueFF <-(1.-(pnorm(-tmp1) + exp(tmp3) *pnorm(-tmp2)))*exp(-k*xx)   #1-G
	}
	if ( all(is.infinite(valueFF)) ) {
		warning(message="Inelegant exit caused by overflow in evaluation of survival function.  
		Check for right-censored data. Try other initial values.")
	}
	
	return(valueFF)	
}

#=========================================================================================================

"survProbInc"<-
function(r,s,k,xx1,xx2)      # calculates incremental survival probability
{
	value.iSP _-(SurvFn(xx2,r,s,k) - SurvFn(xx1,r,s,k))
	value.iSP[value.iSP < 1e-18] _1e-18   # safeguards against taking Log(0)
	value.iSP
}

#=========================================================================================================

"logLikelihood"<-
function(r,s,k,xx1,xx2,Ni)
{
	#returns vector of terms in log likelihood (sum them to get the log likelihood)
	
	# --calculate incremental survival probability--- (safeguraded >1e-18 to prevent log(0))
	iSP <-  survProbInc(r,s,k,xx1,xx2)
		
   loglklhd _ -Ni*log(iSP)
	# add smooth penalty to log-likelihood if k <0
	if (k < 0) {
		loglklhd _loglklhd + k*k*1e4
	}
   return(loglklhd)
}

#=========================================================================================================

"indexFinder"<-
function(x,val)  
#returns the index of the first value of vector x that is <= val.
# returns -1 if no value satisfies the condition
{                 
	idx _(1:length(x))[x<=val][1]
	if(is.na(idx)) {idx_-1}
	idx
}

#=========================================================================================================
"stdErr"<-
function(r,s,k,x1,x2,Ni,pop)
{	
	# function to compute standard error for MLE parameters r,s,k in the vitality model
	# Arguments:
	#  r,s,k - final values of MLE parameters
	#  x1,x2  time vectors (steps 1:(T-1) and 2:T)
	#  Ni - survival fraction
	#  pop - total population of the study
	#
	# Return:
	#  standard error for r,s,k
	# Note: if k <or= 0, can not find std Err for k.
	#
	LL <-function(a,b,c,r,s,k,x1,x2,Ni) {sum(logLikelihood(r+a,s+b,k+c,x1,x2,Ni))}
	
	#initialize hessian for storage
	if (k > 0) {
		hess _matrix(0,nrow=3,ncol=3)
	}
	else {
		hess _matrix(0,nrow=2,ncol=2)
	} 
	
	#set finite difference intervals
	h _.001
	hr _abs(h*r)
	hs _h*s*.1
	hk _h*k*.1


	
	#Compute second derivitives (using 5 point)
	# LLrr
	f0 _LL(-2*hr,0,0,r,s,k,x1,x2,Ni)
	f1 _LL(-hr,0,0,r,s,k,x1,x2,Ni)
	f2 _LL(0,0,0,r,s,k,x1,x2,Ni)
	f3 _LL(hr,0,0,r,s,k,x1,x2,Ni)
	f4 _LL(2*hr,0,0,r,s,k,x1,x2,Ni)

	fp0 _(-25*f0 +48*f1 -36*f2 +16*f3 -3*f4)/(12*hr)
	fp1 _(-3*f0 -10*f1 +18*f2 -6*f3 +f4)/(12*hr)
	fp3 _(-f0 +6*f1 -18*f2 +10*f3 +3*f4)/(12*hr)
	fp4 _(3*f0 -16*f1 +36*f2 -48*f3 +25*f4)/(12*hr)
	
	LLrr _(fp0 -8*fp1 +8*fp3 -fp4)/(12*hr)

	# LLss
	f0 _LL(0,-2*hs,0,r,s,k,x1,x2,Ni)
	f1 _LL(0,-hs,0,r,s,k,x1,x2,Ni)
	# f2 as above
	f3 _LL(0,hs,0,r,s,k,x1,x2,Ni)
	f4 _LL(0,2*hs,0,r,s,k,x1,x2,Ni)

	fp0 _(-25*f0 +48*f1 -36*f2 +16*f3 -3*f4)/(12*hs)
	fp1 _(-3*f0 -10*f1 +18*f2 -6*f3 +f4)/(12*hs)
	fp3 _(-f0 +6*f1 -18*f2 +10*f3 +3*f4)/(12*hs)
	fp4 _(3*f0 -16*f1 +36*f2 -48*f3 +25*f4)/(12*hs)
	
	LLss _(fp0 -8*fp1 +8*fp3 -fp4)/(12*hs)
	
	# LLkk
	if (k > 0) {
		f0 _LL(0,0,-2*hk,r,s,k,x1,x2,Ni)
		f1 _LL(0,0,-hk,r,s,k,x1,x2,Ni)
		# f2 as above
		f3 _LL(0,0,hk,r,s,k,x1,x2,Ni)
		f4 _LL(0,0,2*hk,r,s,k,x1,x2,Ni)

		fp0 _(-25*f0 +48*f1 -36*f2 +16*f3 -3*f4)/(12*hk)
		fp1 _(-3*f0 -10*f1 +18*f2 -6*f3 +f4)/(12*hk)
		fp3 _(-f0 +6*f1 -18*f2 +10*f3 +3*f4)/(12*hk)
		fp4 _(3*f0 -16*f1 +36*f2 -48*f3 +25*f4)/(12*hk)
	
		LLkk _(fp0 -8*fp1 +8*fp3 -fp4)/(12*hk)
	}
	
	
	#-------end second derivs---
	# do mixed partials (4 points)
	# LLrs
	m1 _LL(hr,hs,0,r,s,k,x1,x2,Ni)
	m2 _LL(-hr,hs,0,r,s,k,x1,x2,Ni)
	m3 _LL(-hr,-hs,0,r,s,k,x1,x2,Ni)
	m4 _LL(hr,-hs,0,r,s,k,x1,x2,Ni)

	LLrs _(m1 -m2 +m3 -m4)/(4*hr*hs)
	
	if (k > 0) {
		# LLrk
		m1 _LL(hr,0,hk,r,s,k,x1,x2,Ni)
		m2 _LL(-hr,0,hk,r,s,k,x1,x2,Ni)
		m3 _LL(-hr,0,-hk,r,s,k,x1,x2,Ni)
		m4 _LL(hr,0,-hk,r,s,k,x1,x2,Ni)

		LLrk _(m1 -m2 +m3 -m4)/(4*hr*hk)
	
		# LLsk
		m1 _LL(0,hs,hk,r,s,k,x1,x2,Ni)
		m2 _LL(0,-hs,hk,r,s,k,x1,x2,Ni)
		m3 _LL(0,-hs,-hk,r,s,k,x1,x2,Ni)
		m4 _LL(0,hs,-hk,r,s,k,x1,x2,Ni)

		LLsk _(m1 -m2 +m3 -m4)/(4*hs*hk)
	}

	if (k > 0) {
		diag(hess) _c(LLrr,LLss,LLkk)*pop
		hess[2,1]_hess[1,2]_LLrs*pop
		hess[3,1]_hess[1,3]_LLrk*pop
		hess[3,2]_hess[2,3]_LLrk*pop
	}
	else {
		diag(hess) _c(LLrr,LLss)*pop
		hess[2,1]_hess[1,2]_LLrs*pop
	}
	#print(hess)
	hessInv _solve(hess)
	#print(hessInv)
	
	#compute correlation matrix:
	sz _3
	if (k <= 0) { sz _2 }
	corr_ matrix(0,nrow=sz,ncol=sz)
	for (i in 1:sz) {
		for (j in 1:sz) {
			corr[i,j] _ hessInv[i,j]/sqrt(abs(hessInv[i,i]*hessInv[j,j]))
		}
	}
	#print(corr)
	if ( abs(corr[2,1]) > .98 ) {
		warning("
		WARNING: parameters r and s appear to be closely correlated for this data set.  
		   s.e. may fail for these parameters.     ")
	}
	if (  sz == 3 && abs(corr[3,2]) > .98 ) {
		warning("
		WARNING: parameters s and k appear to be closely correlated for this data set.  
		   s.e. may fail for these parameters.     ")
	}
	if (  sz == 3 && abs(corr[3,1]) > .98 ) {
		warning("
		WARNING: parameters r and k appear to be closely correlated for this data set.  
		   s.e. may fail for these parameters.     ")
	}

	
	se _sqrt(diag(hessInv))
	
	#  Approximate s.e. for cases where calculation of s.e. failed:
	if( sum( seNA_is.na(se) ) > 0 ) {
		se12 _sqrt(diag(solve(hess[c(1,2)	,c(1,2) ])))
		if (k > 0) {
			se13 _sqrt(diag(solve(hess[c(1,3)	,c(1,3) ])))
			se23 _sqrt(diag(solve(hess[c(2,3)	,c(2,3) ])))
		}
		
		if(seNA[1] == T) {
			if(!is.na(se[1]_se12[1]) || (k>0 && !is.na(se[1]_se13[1])) )
				warning("  * s.e. for parameter r is approximate.    ")
			else warning("  * unable to calculate or approximate s.e. for parameter r.    ")
		}
		if(seNA[2] == T) {
			if(!is.na(se[2]_se12[2]) || (k>0 && !is.na(se[2]_se23[1])) )
				warning("  * s.e. for parameter s is approximate.    ")
			else warning("  * unable to calculate or approximate s.e. for parameter s.    ")
		}
		if(k>0 && seNA[3] == T) {
			if(!is.na(se[3]_se13[2]) ||  !is.na(se[3]_se23[2]) )
				warning("  * s.e. for parameter k is approximate.    ")
			else warning("  * unable to calculate or approximate s.e. for parameter k.    ")
		}

	
	
	
	}
		

	#######################
	if (k <= 0) {
		se _c(se,NA)
	}
	return(se)
}
	


	
#=========================================================================================================

"C1.calc"<-
function(pop, isp, Ni)
{
	#  Routine to calculate goodness of fit (Pearson's C -type test)
	#  pop - population number of sample
	#  isp  - Modeled: incemental survivor Probability 
	#  Ni   - Data:  incemental survivor fraction (prob)
	#
	# Returns:  a list containing:
	# C1, dof, Chi2    (retrieve each from list as ..$C1   etc.)
	
if (pop < 35) {
	if (pop < 25) {
		warning(paste("WARNING: sample population (",as.character(pop),") is too small for 
		meaningful goodness of fit measure.  Goodness of fit not being computed"))
		return()
	}
	else {
	warning(paste("WARNING: sample population (",as.character(pop),") may be too small for 
	meaningful goodness of fit measure"))
   }
}
	np <- pop * isp # modeled population at each survival probability level.
	tmpC1 <- 0
	i1 _1; i_1; cnt_0
	len <- length(np)
	while(i <= len) {
		idx <- i1:i
		# It is recommended that each np[i] >5 for meningful results.  Where np<5
		#   points are grouped to attain that level.   I have fudged it to 4.5 ... 
		#  (as some leeway is allowed, and exact populations are sometimes unknown).
		if(sum(np[idx]) > 4.5) {
			cnt _cnt+1
			# Check if enough points remain.  If not, they are glommed onto previous grouping.
			if(i < len && sum(np[(i + 1):len]) < 4.5) {
				idx <- i1:len
				i <- len
			}

			sNi <- sum(Ni[idx])
			sisp <- sum(isp[idx])
			tmpC1 <- tmpC1 + (pop * (sNi - sisp)^2)/sisp
			i1 _i+1			
		}
		
		i _i+1
	}
	C1 <- tmpC1
	dof _cnt-1-3    # degrees of freedom (3 is number of parameters).
		if (dof < 1) {
		warning(paste("WARNING: sample population (",as.character(pop),") is too small for 
		meaningful goodness of fit measure (DoF<1).  Goodness of fit not being computed"))
	   return()
   }

	chi2_qchisq(.95,dof)
	return(C1,dof,chi2)
}

#=========================================================================================================

"plotting"<-
function(r.final,s.final,k.final,mlv,time,sfract,x1,x2,Ni,pplot,tlab,lplot,cplot,gfit,rc.data)
{
	#  Function to provide plotting and goodness of fit computations
	#
	
# --plot cumulative survival---
	if (pplot != F) {
		ext_max(pplot,1)
   		par(mfrow=c(1,1))
   		len_length(time)
   		tmax _ext*time[len]
   		plot(time,sfract,xlab=tlab,ylab="survival fraction",ylim=c(0,1),xlim=c(0,tmax))
		xxx_seq(0,tmax,length=200)
   		lines(xxx,SurvFn(xxx,r.final,s.final,k.final))
   		title("Cumulative Survival Data and Vitality Model Fitting")

       # --likelihood and likelihood contour plots---
		if(lplot != F) {
			profilePlot(r.final,s.final,k.final,x1,x2,Ni,mlv,cplot)
		}
	}

# --calculations for goodness of fit---
	if(gfit!=F) { # then gfit must supply the population number
   		isp _survProbInc(r.final,s.final,k.final,x1,x2)
		C1dof <- C1.calc(gfit,isp,Ni)
		C1 _C1dof$C1
		dof _C1dof$dof
		chi2 _C1dof$chi2
				
		print(paste("Pearson's C1=",as.character(round(C1,3)),"   chisquared =",
				as.character(round(chi2,3)),"on",as.character(dof),"degrees of freedom "))
		
		#   Note: The hypothesis being tested is whether the data could reasonably have come 
		#   from the assumed (vitality) model.
		if (C1 > chi2) {
			print("C1 > chiSquared;   should reject the hypothesis becasue C1 falls outside the 95% confidence interval.")
		}
		else {
			print("C1 < chiSquared;   should Not reject the hypothesis because C1 falls inside the 95% confidence interval")		
		}
	}
		
	# --Incremental mortality plot
	if (pplot != F) {
		par(mfrow=c(1,1))
		if (rc.data != F) {
			ln _length(Ni)-1
			x1 _x1[1:ln]
			x2 _x2[1:ln]
			Ni _Ni[1:ln]
		}
		ln _length(Ni)
		#scale_(x2-x1)[Ni==max(Ni)]
		scale_max( (x2-x1)[Ni==max(Ni)] )
					
		ext_max(pplot,1)
		
		npt_200*ext
		xxx _seq(x1[1],x2[ln]*ext,length=npt)
		xx1 _xxx[1:(npt-1)]
		xx2 _xxx[2:npt]
		sProbI _survProbInc(r.final[1],s.final[1],k.final[1],xx1,xx2)
			
		ytop _1.1*max( max(sProbI/(xx2-xx1)),Ni/(x2-x1) )*scale		
		plot((x1+x2)/2,Ni*scale/(x2-x1),ylim=c(0,ytop),xlim=c(0,ext*x2[ln]),xlab=tlab,ylab="incremental mortality")
		title("Probability Density Function")
		lines((xx1+xx2)/2,sProbI*scale/(xx2-xx1))
	}
	
	return()	 
}

#=========================================================================================================

"profilePlot"<-
function(r.f,s.f,k.f,x1,x2,Ni,mlv,cplot) 
{
	 #
	 # mlv = value of max likelihood 
	
    

# likelihood plots
	# These plots are really not "profile" plots in that while one parameter is varied, 
	# the others are held fixed, not re-optimized. 
	
	SLL <- function(r,s,k,x1,x2,Ni) { sum(logLikelihood(r,s,k,x1,x2,Ni)) }

	rf_.2; sf_.5; kf_1.0; fp_40 # rf,sf,kf - set profile plot range (.2 => plot +-20%), 2*fp+1 points
	rseq _seq((1-rf)*r.f,(1+rf)*r.f, (rf/fp)*r.f)
	sseq _seq((1-sf)*s.f,(1+sf)*s.f, (sf/fp)*s.f)
	if (k.f > 0) {
		kseq _seq((1-kf)*k.f,(1+kf)*k.f, (kf/fp)*k.f)
	}
	else {  #if k=0..
		kseq _seq(.00000001,.1,length=(2*fp+1))
	}

	rl _length(rseq)
	tmpLLr _rep(0,rl)
	tmpLLs _tmpLLr
	tmpLLk _tmpLLr
	for (i in 1:rl) {
		tmpLLr[i] _SLL(rseq[i],s.f,k.f,x1,x2,Ni)
		tmpLLs[i] _SLL(r.f,sseq[i],k.f,x1,x2,Ni)	
		tmpLLk[i] _SLL(r.f,s.f,kseq[i],x1,x2,Ni)
	}
	
	par(mfrow=c(2,3))
	
	rlim1 _rseq[1]
	rlim2 _rseq[rl]
	if (r.f < 0) {  #even though r should not be <0
		rlim2 _rseq[1]
		rlim1 _rseq[rl]
	}
			  
	plot(r.f,LL_SLL(r.f,s.f,k.f,x1,x2,Ni),
			xlim=c(rlim1,rlim2), xlab="r",ylab="Likelihood");  
	lines(rseq,tmpLLr)
	plot(s.f,LL,xlim=c(sseq[1],sseq[rl]), xlab="s",ylab="Likelihood");  
	lines(sseq,tmpLLs)
	title("Likelihood Plots")
	plot(k.f,LL,xlim=c(kseq[1],kseq[rl]), 
	       ylim=c(1.1*min(tmpLLk)-.1*(mLk_max(tmpLLk)),mLk),xlab="k",ylab="Likelihood"); 
	lines(kseq,tmpLLk)
			
			
	# --  for contour plotting  ----------------------------------------
		
 	if (cplot==T) {
		rl2_(rl+1)/2; rl4_20; st_rl2-rl4; nr_2*rl4+1
		tmpLLrs _matrix(rep(0,nr*nr),nrow_nr,ncol=nr)
		for(i in 1:nr) {
			for(j in 1:nr) {
				tmpLLrs[i,j] _SLL(rseq[i+st-1],sseq[j+st-1],k.f,x1,x2,Ni)
			}
		}
			
		lvv_seq(mlv,1.02*mlv,length=11)  #99.8%, 99.6% ... 98%
		par(mfrow=c(1,1))
		contour(rseq[st:(rl2+rl4)],sseq[st:(rl2+rl4)],tmpLLrs,levels=lvv,xlab="r",ylab="s")
		title("Likelihood Contour Plot of r and s.",  "Outermost ring is likelihhod 98% (of max) level, innermost is 99.8% level.")
		points(r.f,s.f,pch="*",cex=3.0)
		points(c(rseq[st],r.f),c(s.f,sseq[st]),pch="+",cex=1.5)
	}

}

#=========================================================================================================