.packageName <- "accuracy"
# frexp
#
# Extract mantissa's and exponents from a vector of numbers
# R wrapper around the frexp() C-library call.
#
# Part of the Accuracy package. Available from www.r-project.org and
# www.hmdc.harvard.edu/numerical_issues/
#
#    Copyright (C) 2004  Micah Altman
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


"frexp" <-
function(v) {
	r = replicate(length(v), NA)
	r = cbind(r,r)
	i = which(is.finite(v) | !is.na(v))
	tmp = .C("R_frexp", PACKAGE="accuracy", as.double(v[i]), as.integer(length(v[i])),
		mantissa=double(length(v[i])),
		exponent=integer(length(v[i])))
	r[i,] = cbind(tmp$mantissa,tmp$exponent)
	dimnames(r)[[2]]=c("Mantissa","Exponent")
	return(r)
}


"frexpTest" <-
function(silent=TRUE) {
	d=options()$digits
	options(digits=15)
	f=round(frexp(c(1,2,4,1.1,2.1,4.1,1000,20000,1.02E213)),digits=5)
	x=cbind(rbind(.5,.5,.5,.55,.525,.5125,.97656,.61035,.75747),
		 rbind(1,2,3,1,2,3,10,15,708))
	options(digits=d)
	ret=(sum(f==x)==18)
	if (!ret && !silent) {
		warning("Failed frexp self test.")
	}
	return(ret)
}

"LRE"<-
function(x,correct, use.LAE=TRUE) {
        zeros = which(x==0)
        nonzeros = which(x!=0)
        nas = which(x==NA)

        res=double(length(x))

        res[nas]=NA
        res[nonzeros]= -1 * 
                log10(abs(x[nonzeros]-correct[nonzeros])/correct[nonzeros]) 
        if (use.LAE) {
                res[zeros]=-1*log10(abs(x[zeros]-correct[zeros]))
        } else {
                res[zeros]=NA
        }

        return(res)
}

# globaltests
#
# Tests for global optimality of non-linear and maximum-likelihood solutions
#
# Part of the Accuracy package. Available from www.r-project.org and
# www.hmdc.harvard.edu/numerical_issues/
#
#    Copyright (C) 2004  Micah Altman , Michael McDonald
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


dehaan<-function(llTest, llMax, pval=.05 ) {
  # returns TRUE if "llMax" likelihood is greater than
  # the (1-pval) confidence interval for the true optimum
  # derived from "llTest" likelihoods.

  if ( !is.numeric(llTest) ||  !is.numeric(llMax) ) {
        stop("llTest must be numeric")
  }

  if (length(llTest)<2) {
        stop("Size of llTest must be greater than 1")
  }

  if (length(pval) !=1 || pval>=1 || pval<0)  {
        stop("pval must be scalar: 0< pval <1")
  }

  x = sort(llTest)
  n = length(llTest)


  # De Hann (1981) proves this test statistic works for some k(n) where
  # k(n)/n->0 as n->infinity.  A likely candidate is k=sqrt(n)

  k = n^(1/2)
  alpha = k/2

  lp = x[n]+ ( (x[n])-(x[n-1]) ) / (pval ^ (-1/alpha) -1)

  if (lp>llMax)  {
    return(FALSE)
  } else  {
    return(TRUE)
  }
}

dehaanSelftest<-function(silent = TRUE) {
        
        # Tests of deHaan function

        ret = TRUE

        # This test data should always return a FALSE.

        l1= c(1,2,3,4,5)
        max1 = 4; 

        if (dehaan(c(1,2,3,4,5),4) ||
	    !dehaan(c(1,2,3,4,5),1000) )  {
                ret = FALSE
                 if (!silent) {
                        warning("failed selftest ")
                }
        }

	if (  dehaan(1:100,1) ||
	      dehaan(1:100,101) || 
	      dehaan(1:100,99) || 
	      dehaan(1:20,21,pval=.9999) ||
              !dehaan(1:100,101,pval=.001)|| 
	      !dehaan(1:20,21) 
 	) {
           ret = FALSE
           if (!silent) {
                  warning("failed selftest")
           }
        } 

        return(ret)
}


starr<-function(betas, tol=.0001, dmethod="euclidean") {

  if (!is.matrix(betas)) {
    stop("betas must be matrix")
  }
  if (length(tol)!=1 || (tol<0) || !is.numeric(tol)) {
    stop("tolmust be scalar >=0")
  }

        return(ret)
}


starr<-function(betas, tol=.0001, dmethod="euclidean") {

  if (!is.matrix(betas)) {
    stop("betas must be matrix")
  }
  if (length(tol)!=1 || (tol<0) || !is.numeric(tol)) {
    stop("tolmust be scalar >=0")
  }

  # find unique (within tolerance) optima, count them

  n = nrow(betas)
  if (n==1)  {
	return(1)
  }
 
  dm = as.matrix(dist(betas,method=dmethod))
  optc = integer(n)
  for (i in 1:n) {
        optc[i] = sum(dm[i,]<tol)
  }

  sortInd = sort(optc,decreasing=TRUE,index.return=TRUE)$ix
  for (i in sortInd) {
    if (optc[i]>0) {
    tmp=optc[i]
        optc[dm[i,]<tol]=0; 
        optc[i] = tmp
    }
  }
  nopt = sum(optc>0); 

  ndouble = sum(optc==2)
  nsingle = sum(optc==1)

  rv = nsingle/n + 2*ndouble/(n*(n-1))
  
  return(rv)
                              
}

starrSelftest<-function (silent=TRUE) {

  # BOD test
  #
  # For BOD example, see help("BOD")

  ret = TRUE
  

  x=rbind(c(1,1,1), c(1,2,1), c(1,1.1,1), c(1,2,1), c(3,4,5))
  if ( (starr(x)!=.7) || 
       (starr(rbind(1)) !=1) ||
       (starr(rbind(1,1,2,2)) !=1/3) ||
       (starr(rbind(1,1.00001,2,2)) !=1/3) ||
       (starr (rbind(1,1,1,1,1,1,1,1,1,1,1,1)) !=0) 
     ) 
   {
    	ret=FALSE
    	if (!silent) {
       		warning("Starr test for global optimum failed selftest")
    	}
   }
  return(ret)

}

starrRun<-function(start,optfunc,...) {
  r= NULL
  dstart = as.data.frame(start)
  for (i in 1:nrow(dstart)) {
    r = rbind(r, coef(optfunc(..., t(dstart)[,i])))
  }
}
# perturb.R
#
# Performs sensitivity analysis of non-linear and linear models
# through data pertubations
#
# Part of the Accuracy package. Available from www.r-project.org and
# www.hmdc.harvard.edu/numerical_issues/
#
#    Copyright (C) 2004  Micah Altman
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



# perturb
#
# Replicates an analysis with randomly perturbed data
#
# data = data
# statistic = stat function to run
# ptb.R = number of replications
# ran.gen= function or vector of functions to apply to data to add noise
# ptb.s = size, or vector of sizes for noise generators

perturb<-function(data,statistic,..., ptb.R=1,ptb.ran.gen=PTBms,
	ptb.s=NULL) {
	ptb.R=as.integer(trunc(ptb.R))

	if ( (is.integer(ptb.R)==FALSE) || (ptb.R<=0) || (length(ptb.R)>1) ) {
		stop("ptb.R must be int > 0")
	}
	retval = list(ptb.R);
	for (i in 1:ptb.R) {
		if (is.null(ptb.s)) {
		   retval[[i]]=
		    perturbHarness(data=data,ran.gen= ptb.ran.gen,statistic=statistic,...)
		} else  {
		   retval[[i]]=
		    perturbHarness(data=data,ran.gen= ptb.ran.gen,statistic=statistic,...,
			ptb.s=ptb.s)
		} 
	}

	# replicate swallows ... in R 1.9.0 ???!!??, worked in 1.6-1.8
	#retval = replicate(ptb.R,
	#	perturbHarness(data=data,ran.gen= ptb.ran.gen,statistic=statistic,...)
	#, simplify=FALSE);

	attr(retval,"ran.gen")= ptb.ran.gen
	attr(retval, "R") = ptb.R
	attr(retval, "s") = ptb.s
	attr(retval, "statistic") = statistic
	class(retval)="perturb"
	return(retval)
}

print.perturb<-function(x,...) {
	cat("Replications: \n",attr(x,"R"),"\n\n")
	cat("ran.gen: \n")
	print(attr(x,"ran.gen"))
	cat("s: \n")
	print(attr(x,"s"))
	cat("statistic: \n")
	print(attr(x,"statistic"))
}

perturbHarness<-function(data,ran.gen,statistic,..., ptb.s=NULL) {
	ndata = as.data.frame(data)
	if (length(ran.gen)==1) {
		ind=which(sapply(ndata,is.numeric))
		ran.gen=replicate(ncol(ndata),ran.gen)
		if (!is.null(ptb.s)) {
			ptb.s = replicate(ncol(ndata),ptb.s)
		}
	} else if (length(ran.gen)!=ncol(ndata)) {
		stop("ran.gen must be a single function, or a vector of functions of length ncol(df)")
	} else {
		ind = which(sapply(ran.gen,is.function))
	}

	for ( i in ind ) {
		if (is.null(ptb.s)) {
			ndata[i] = ran.gen[[i]](ndata[i])
		} else {
			ndata[i] = ran.gen[[i]](ndata[i],size=ptb.s[i])
		}
	}
	
	stat = statistic(data=ndata,...)
	
	return(stat)
}


############################################################
# 
# Perturbations for vectors
#
# These functions take a vector and apply a mean-zero
# random perturbation to them.
#
# Fairly forgiving as to inputs. Will accept a list,
# matrix, scalar, or dataframe as well. However,
# applying a centered perturbation to a dataframe or
# matrix centers it on the whole dataframe or matrix, 
# rather than centering it on the vector. This is probably
# not what you want.
#
#
# Perturb.unif -- uniform perturbations 
# Perturb.meps -- perturbations on the order of storage
#		  roundoff error
# Perturb.norm -- normal perturbations
# 
#
# Inputs: 
#	x - vector or matrix
#	s - size of perturbation
#	centered - boolean, center disturbance to guarantee
#		(up to numeric tolerance) mean 0 for uniform
#		and normal disturbances, when used by itself
#	scaled - make disturbance size scaled to observation
#		size
#	Note: Using both scaled=TRUE and centered=TRUE, 
#	 	cannot guarantee mean 0-- the centering
# 		occurs before the scaling of the noise.
#		Once rescaled, the noise is no longer
#		guaranteed to be mean zero.

#
# Returns:
#	 perturbed vector, matrix or dataframe
#
##############################################################

PTBunif<-function(x ,size=1 , centered=FALSE, scaled=FALSE) {
	delta=runif(length(as.matrix(x)))*size; 
	
	if (scaled && centered) {
		warning("Using scaled and centering together is rarely what you really want to do.")
	}

	# centering to guarantee mean 0
	if (centered) {		
		delta=delta-mean(delta)
	}

	# scaled disturbance
	if (scaled) {
		delta = delta * x
	}	
	
	return(x+delta)
}

PTBnorm<-function(x , size=1, centered=FALSE, scaled=FALSE) {
	delta = rnorm( length(as.matrix(x)), mean=0, sd= 1) * size
	

	if (scaled && centered) {
		warning("Using scaled and centering together is rarely what you really want to do.")
	}

	# centering to guarantee mean 0
	if (centered) {
		delta = (delta - mean(delta))/sqrt(var(delta))
	}

	if (scaled) {
		delta = delta * x
	}	
	
	return(x+delta)
}

PTBmeps<-function(x, centered=FALSE, scaled=FALSE, size=1) {	
	# scaled=FALSE leads to larger roundoff 
	if (!scaled) {
		warning("scaled=FALSE probably not what you want for PTBmeps")
	}

	n = length(as.matrix(x))
	
	# note centering here does not guarantee zero mean!
	if (centered) {
		delta=integer(n)
		delta[1:floor((n+1)/2)]=.Machine$double.eps *size
		delta[ceiling((n+1)/2):n]=-1*.Machine$double.neg.eps *size
		delta=sample(delta,n)
	} else {
		delta= runif(n)*2-1
		delta= floor(delta)*.Machine$double.neg.eps +
			ceiling(delta)*.Machine$double.eps
	}

	if (scaled) {
		delta = delta*x
	}
	return(x+delta)
}

# Wrapper functions for bounded perturbations

PTBmsb<-function(x,size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="trunc",scaled=TRUE,
		ran.gen=PTBmeps)
}

PTBmsbr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="resample", scaled=TRUE,
		ran.gen=PTBmeps)
}

PTBubr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="resample",
		ran.gen=PTBunif)
}

PTBubrr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="relresample",
		ran.gen=PTBunif)
}

PTBnbr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="resample",
		ran.gen=PTBunif)
}

PTBnbrr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="relresample",
		ran.gen=PTBunif)
}

PTBusbr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="resample",scaled=TRUE,
		ran.gen=PTBunif)
}

PTBusbrr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="relresample",scaled=TRUE,
		ran.gen=PTBunif)
}

PTBnsbr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="resample", scaled=TRUE,
		ran.gen=PTBunif)
}

PTBnsbrr<-function(x, size=1, lbound=0,ubound=1) {
	ptbBndHarness(x,lbound,ubound,size=size,mode="relresample",scaled=TRUE,
		ran.gen=PTBunif)
}

# bounded perturbation -- takes set of samples and ensure results
# are within variable bounds

ptbBndHarness<-function(x,lbound,ubound,s=NULL,mode="trunc",
	ran.gen=PTBunif, maxiter=5000,...) {

	# check args
	if (mode!="trunc" && mode !="resample" && mode != "relresample") {
		stop("unknown mode")
	}

	# initialize bounds vectors
	if (length(lbound)==1) {
		lbound=replicate(length(x),lbound)
	} else if (length(lbound)!=length(x)) {
		stop("lbound must be scalar or same length as x")
	}
	if (length(ubound)==1) {
		ubound=replicate(length(x),ubound)
	} else if (length(ubound)!=length(x)) {
		stop("ubound must be scalar or same length as x")
	}

	# set sizes for relresample mode
	if (mode=="relresample") {
		tmps=pmin(x-lbound,ubound-x)
		s = pmin(s,tmps)
	}


	# initial sample	
	if (is.null(s)) {
		retv = ran.gen(x,...)
	} else {
		retv = ran.gen(x, s=s,...)
	}
	
	if ( mode=="resample" || mode=="relresample") {
		badi = which( retv<lbound | retv>ubound )
		iter=0
		while (length(badi)>0 && iter<=maxiter) {
			if (is.null(s)) {
				retv[badi] = ran.gen(x[badi],...)
			} else {
				retv[badi] = ran.gen(x[badi], s=s,...)
			}
			badi = which( retv<lbound | retv>ubound )
			iter=iter+1;	
		}
		if (iter>maxiter) {
			warning("resampling: maximum iterations exceeded, truncating to bounds")
		}
	} 

	# this is 'trunc' mode, also catches exceeded iterations on
	# resampling, and numerical issues in comparisons

	retv=pmax(lbound,retv)
	retv=pmin(ubound,retv)

	return(retv)
	
}


###############################################################
#
# These are wrappers around ptb.{meps,unif,norm}
# to make it easier to use the ptb vector
# functions from the perturb() harness.
#
###############################################################

PTBi<-function(x,size=1) {
	return(x)
}

PTBus<-function(x ,size=1) {
	return(PTBunif(x,size=size,scaled=TRUE))
}

PTBuc<-function(x ,size=1) {
	return(PTBunif(x,centered=TRUE,size=size))
}

PTBu<-function(x ,size=1) {
	return(PTBunif(x,centered=FALSE,scaled=FALSE,size=size))
}


PTBns<-function(x ,size=1) {
	return(PTBnorm(x,scaled=TRUE,size=size))
}

PTBnc<-function(x ,size=1) {
	return(PTBnorm(x,centered=TRUE,size=size))
}

PTBn<-function(x ,size=1) {
	return(PTBnorm(x,centered=FALSE,scaled=FALSE,size=1))
}

PTBms<-function(x,size=1) {
	return(PTBmeps(x,centered=FALSE,scaled=TRUE,size=size))
}


#
# Summary Functions
#

summary.perturb<-function (object,...) {
	n = length(object)
	s = vector(n,mode="list")

	# generate individual summaries for list of replications
	for (i in 1:n) {
		tmp = summary(object[[i]])
		if (inherits(tmp,"summary.lm")) {
			coef.names = dimnames(tmp[["coefficients"]])[[1]]
			coef.betas = tmp[["coefficients"]][,1]
			coef.stderrs = tmp[["coefficients"]][,2]
			coef.formula = tmp[["call"]]
		} else if (inherits(tmp,"summary.mle")) {
			coef.names = dimnames(attr(tmp,"coef"))[[1]]
			coef.betas = attr(tmp,"coef")[,1]
			coef.stderrs = attr(tmp,"coef")[,2]
			coef.formula = attr(tmp,"call")
		} else if (inherits(tmp,"summary.nls")) {
			coef.names = dimnames(tmp[["parameters"]])[[1]]
			coef.betas = tmp[["parameters"]][,1]
			coef.stderrs = tmp[["parameters"]][,2]
			coef.formula = tmp[["formula"]]
		} else if (inherits(tmp,"summary.glm")) {
			coef.names = dimnames(tmp[["coefficients"]])[[1]]
			coef.betas = tmp[["coefficients"]][,1]
			coef.stderrs = tmp[["coefficients"]][,2]
			coef.formula = tmp[["terms"]]
		} else {
			coef.names=names(coef(tmp))
			coef.betas=coef(tmp)
			coef.stderrs=NULL
			coef.formula=NULL
			if (is.null(coef.betas)) {
				stop("Don't know how to summarize replications of type ",
					class(tmp) )
			}
		}
	
 		attr(tmp,"coef.names") = coef.names
		attr(tmp,"coef.betas") = coef.betas
		attr(tmp,"coef.stderrs") = coef.stderrs
		attr(tmp,"coef.formula") = coef.formula
		s[[i]]=tmp
	}
	
	# check consistency and summarize
	coef.names = attr(s[[1]],"coef.names")
	coef.names.m = coef.names
	coef.betas = attr(s[[1]],"coef.betas")
	coef.betas.m = coef.betas
	coef.stderrs= attr(s[[1]],"coef.stderrs")
	coef.stderrs.m = coef.stderrs;	
	coef.formula = attr(s[[1]],"coef.formula")
	if (n>1) {
	   for (i in 2:n) {
		if ( sum( attr(s[[i]],"coef.names") != coef.names)>0
		   || sum( attr(s[[i]],"coef.formula") != coef.formula)>0
		   || sum( length(attr(s[[i]],"coef.betas")) != length( coef.betas))>0
		   || sum( length(attr(s[[i]],"coef.stderrs")) != length( coef.stderrs))>0
		) {
			warning("replications does not match (", i,")")
		}
		coef.names.m = rbind(coef.names.m, attr(s[[i]],"coef.names"))
		coef.betas.m = rbind(coef.betas.m, attr(s[[i]],"coef.betas"))
		coef.stderrs.m = rbind(coef.stderrs.m, attr(s[[i]],"coef.stderrs"))
	   }
	}


	row.names(coef.betas.m) = NULL
	row.names(coef.stderrs.m) = NULL
	attr(s,"coef.names.m") = coef.names.m
	attr(s,"coef.betas.m") = coef.betas.m
	attr(s,"coef.stderrs.m") = coef.stderrs.m
	attr(s,"coef.formula") = coef.formula


	attr(s, "ran.gen")= attr(object, "ran.gen")
	attr(s, "R") = attr(object, "R")
	attr(s, "s") = attr(object, "s")
	attr(s,"statistic") = attr(object, "statistic")

	class(s)="perturbS";	
	return(s)
}

print.perturbS<-function(x,...) {

	cat("Replications: \n",attr(x,"R"),"\n\n")
	cat("ran.gen: \n")
	print(attr(x,"ran.gen"))
	cat("s: \n")
	print(attr(x,"s"))
	cat("statistic: \n")
	print(attr(x,"statistic"))

	cat("formula: \n","\n\n")
	print(attr(x,"coef.formula"))

	cat("betas:\n\n")
	print(summary(attr(x,"coef.betas.m")))
	cat("stderrs:\n\n")
	print(summary(attr(x,"coef.stderrs.m")))
	
}


perturbTest<-function(silent=TRUE) {
	status=TRUE
	if(version$major<2) {
		data(longley, package="base")
	} else {
		data(longley, package="datasets")
	}

	for (i in c(PTBnc,PTBuc)) {
		pl = i(longley)
		if (sum(pl==longley) != 0) {
			status=FALSE
			if (!silent) {
				warning("not perturbing --")
				print(i)
			}
		}
		if (abs(mean(mean(pl-longley)))>.000001) {
			status=FALSE
			if (!silent) {
				warning("perturbations too big--")
				print(i)
			}
		}
	}
	for (i in c(PTBms,PTBn,PTBu)) {
		pl = i(longley)
		if (sum(pl==longley) != 0) {
			status=FALSE
			if (!silent) {
				warning("not perturbing --")
				print(i)
			}
		}
		if (abs(mean(mean(pl-longley)))>1) {
			status=FALSE
			if (!silent) {
				warning("perturbations too big--")
				print(i)
			}
		}
	}
	for (i in c(PTBns,PTBus)) {
		pl = i(longley)
		if (sum(pl==longley) != 0) {
			status=FALSE
			if (!silent) {
				warning("not perturbing --")
				print(i)
			}
		}
		if (abs(mean(mean(pl-longley)))>abs(mean(mean(longley)))) {
			status=FALSE
			if (!silent) {
				warning("perturbations too big--")
				print(i)
			}
		}
	}
	for (i in c(PTBnbr,PTBubr)) {
		t = runif(20)*2-1
		pl = i(t, lbound=-1, ubound=1)
		if (sum(t>=1) > 0 || sum(t<=-1)> 0) {
			status = FALSE
			if (!silent) {
				warning("bounded perturbations not bounded--")
				print(i)
			}
		}
	}

	#data test
	plongley = perturb(longley,lm,Employed~., ptb.R=10,
    	   ptb.ran.gen=c(PTBi, replicate(5,PTBus),PTBi), ptb.s=c(1,replicate(5,.001),1))
	sp=summary(plongley)
	coef= attr(sp,"coef.betas.m")
	if (nrow(coef)!=10 || ncol(coef)!=7) {
		status = FALSE
		if (!silent) {
			warning("perturb framework malfunction")
			print(i)
		}
	}

	return(status)
}
# sechol.R
#
# Schnable-Eskow generalized cholesky.
#
# Part of the Accuracy package. Available from www.r-project.org and
# www.hmdc.harvard.edu/numerical_issues/
#
#    Copyright (C) 2004  Jeff Gill
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

"sechol" <- function(A, tol = .Machine$double.eps, silent= TRUE )  {
    n <- nrow(A)
    L <- matrix(rep(0,n*n),ncol=ncol(A))
    tau <- tol ^(1/3)  # made to match gauss
    gamm <- max(A)
    deltaprev <- 0
    Pprod <- diag(n)
    if (n > 2)  {
      for (k in 1:(n-2))  {
	if( (min(diag(A[(k+1):n,(k+1):n]) - A[k,(k+1):n]^2/A[k,k]) < tau*gamm) 
		&& (min(svd(A[(k+1):n,(k+1):n])$d)) < 0) {
	    dmax <- order(diag(A[k:n,k:n]))[(n-(k-1))]
	    if (A[(k+dmax-1),(k+dmax-1)] > A[k,k])  {
		if (!silent) {
	           print(paste("iteration:",k,"pivot on:",dmax,"with absolute:",(k+dmax-1)))
		}
	        P <- diag(n)
                Ptemp <-  P[k,]; P[k,] <- P[(k+dmax-1),]; P[(k+dmax-1),] = Ptemp
	        A <- P%*%A%*%P
	        L <- P%*%L%*%P
	        Pprod <- P%*%Pprod
	    }
	    g <- rep(0,length=(n-(k-1)))
	    for (i in k:n)  {
		if (i == 1) sum1 <- 0
		else sum1 <- sum(abs(A[i,k:(i-1)]))
		if (i == n) sum2 <- 0
		else sum2 <- sum(abs(A[(i+1):n,i]))
		g[i-(k-1)] <- A[i,i] - sum1 - sum2
	    }
	    gmax <- order(g)[length(g)]
	    if (gmax != k)  {
		if (!silent) {
	            print(paste("iteration:",k,
			"gerschgorin pivot on:",gmax,"with absolute:",(k+gmax-1)))
		}
	        P <- diag(ncol(A))
                Ptemp <-  P[k,]; P[k,] <- P[(k+dmax-1),]; P[(k+dmax-1),] = Ptemp
	        A <- P%*%A%*%P
	        L <- P%*%L%*%P
	        Pprod <- P%*%Pprod
	    }
	    normj <- sum(abs(A[(k+1):n,k]))
	    delta <- max(0,deltaprev,-A[k,k]+max(normj,tau*gamm))
	    if (delta > 0)  {
		A[k,k] <- A[k,k] + delta
		deltaprev <- delta
	    }
	}
        L[k,k] <- A[k,k] <- sqrt(A[k,k])
	for (i in (k+1):n)  {
	    L[i,k] <- A[i,k] <- A[i,k]/L[k,k]
	    A[i,(k+1):i] <- A[i,(k+1):i] - L[i,k]*L[(k+1):i,k]
	    if(A[i,i] < 0) A[i,i] <- 0
	}
      }
    }
    A[(n-1),n] <- A[n,(n-1)]
    eigvals <- eigen(A[(n-1):n,(n-1):n])$values
    delta <- max(0,deltaprev,
        -min(eigvals)+tau*max((1/(1-tau))*(max(eigvals)-min(eigvals)),gamm))
    if (delta > 0)  {
	if (!silent) {
        	print(paste("delta:",delta))
	}
        A[(n-1),(n-1)] <- A[(n-1),(n-1)] + delta
        A[n,n] <- A[n,n] + delta
        deltaprev <- delta
    }
    L[(n-1),(n-1)] <- A[(n-1),(n-1)] <- sqrt(A[(n-1),(n-1)])
    L[n,(n-1)] <- A[n,(n-1)] <- A[n,(n-1)]/L[(n-1),(n-1)]
    L[n,n] <- A[n,n] <- sqrt(A[n,n] - L[n,(n-1)]^2)
    
   r = t(Pprod)%*%t(L)%*%t(Pprod)
   attr(r,"delta")=delta
   return(r)
}

"secholTest"<-function(silent=TRUE) {
     rv = TRUE
     # non singular
     S <- matrix(c(2,0,2.4,0,2,0,2.4,0,3),ncol=3)
     rv = (sum( signif(chol(S),digits=14) == signif(sechol(S),digits=14)) ==9)
     if (!rv && !silent) {
		warning("sechol alters PD matrix")
     }
     S <- matrix(c(2,0,10,0,2,0,10,0,3),ncol=3)
     t =(
	  sum(
	   signif(sechol(S), digits=10) == 
	   matrix(c(1.414213562,0,0,0,1.414234971,0,7.071067812
		,0,0.007781680058), ncol=3)
	  ) == 9 )

     if (!t && !silent) {
		warning("sechol results don't match benchmark")
     }
     rv= rv && t
     return(rv)
}








# truerandom.R
#
# R methods to obtain true random numbers from the /dev/random
# entropy collector.
#
# Part of the Accuracy package. Available from www.r-project.org and
# www.hmdc.harvard.edu/numerical_issues/
#
#    Copyright (C) 2004  Micah Altman
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#
# Uses Linux /dev/random if available to get high-quality true random
# numbers for use in setting the PRNG seed. Uses hotbits to get
# numbers as well.
#


"runifT" <- function(n, min=0, max=1) {
	if (min>=max) {
		stop("Max must be > min")
	}
	tmp = trueRandom(n)
	if (is.null(tmp)) {
		warning("No entropy available, returning pseudo-random numbers")
		r = runif(n, min, max)
	} else { 
		r =(tmp/.Machine$integer.max  + 1) * ((max-min)/2)
	}
	return(r)
}

"trueRandom" <-function (n) {
	if (length(n)>1) {
		size=length(n)
	} else {
		size = n
	}

	if (!exists(".EntropyPool",envir=.GlobalEnv)) {
		pool=refreshPool(silent=TRUE)
		if (is.null(pool)) {
			return(NULL)
		}
	} else {
		pool=get(".EntropyPool",envir=.GlobalEnv)
	}

	tr = integer(size)
	i = 1
	while (i<=size) {
		if (pool$current < 1)  {
			pool=refreshPool(silent=TRUE)
			next
		} 
		tr[i] = pool$pool[pool$current]
		pool$current = pool$current -1
		i= i+1
	}	

	assign(".EntropyPool",pool,envir=.GlobalEnv)
	return(tr)
}

"resetSeed" <-function() {
	s = trueRandom(1)
	if (is.null(s)) {
		warning("No entropy available, using system time as seed.")
		s = as.integer(Sys.time())
	}
	set.seed(s)
}

"initPool"<-function(size=512, hbok=TRUE, devrndok=TRUE, silent=FALSE) {
	entropypool = list()
	entropypool$size=size
	entropypool$current=0
	entropypool$pool=integer(length=size)
	class(entropypool)="EntropyPool"
	
	if (size<1 || size > 10000/.Machine$sizeof.long) {
		warning("size out of range")
		return (NULL)
	}

	w = options("warn")
	options(warn=-1)
	if (devrndok) {
	   tri = integer()
           tr = try({tri=readBin('/dev/random', integer(0), signed=FALSE)},
                silent=TRUE)
           if (inherits(tr, "try-error") || (length(tri) == 0)) {
                devrndok=FALSE
           }   
	}
  
	if (hbok) {
	   tri = integer()
	   hb = try({hburl(bytes= .Machine$sizeof.long)})
           if (is.null(hb) || inherits(hb, "try-error")) {
           	tr = try({tri=readBin(hb, integer(0), signed=FALSE)},
                	silent=TRUE)

                 if (inherits(tr, "try-error") || (length(tri) == 0)) {
                   hbok=FALSE
                 }  else {
	   		try(close(hb), silent=TRUE)
		}
	   } else {
                   hbok=FALSE
	   }
	}
	options(warn=as.integer(w))
	entropypool$hbok=hbok
	entropypool$devrndok=devrndok
	if (!devrndok && !hbok ) {
		if (!silent) {
			warning("initialization failed, no true random sources found")
		}
		return(NULL)
	}
	assign(".EntropyPool",entropypool,envir=.GlobalEnv)
	return(entropypool)
}

"hburl" <-function(bytes=1,fmt="bin") {
	hbstring = paste (
		"http://www.fourmilab.ch/cgi-bin/uncgi/Hotbits?"
		,"nbytes=",bytes,"&fmt=",fmt, sep="")
	return(url(hbstring,open="rb"))
}

"refreshPool"<-function(silent=FALSE) {
	if (!exists(".EntropyPool",envir=.GlobalEnv)) {
		if (is.null(initPool(silent=TRUE))) { 
			if (!silent) {
				warning("Could not create pool")
			}
			return(NULL)
		}
	}
	pool=get(".EntropyPool",envir=.GlobalEnv)
	if (pool$devrndok) {
		con = "/dev/random"
		pool$pool = readBin(con, integer(0), signed=FALSE,n=pool$size)
	} else if (pool$hbok) {
		con = hburl(bytes=pool$size)
		pool$pool = readBin(con, integer(0), signed=FALSE,n=pool$size)
		close(con)
	}
	pool$current=length(pool$pool)
	assign(".EntropyPool",pool,envir=.GlobalEnv)
	return(pool)
}
.onLoad <- function(lib, pkg) {
  if((version$major<1) ||(version$minor<6 && version$major==1)) {
    stop("This version for R 1.6 or later")
  }
  if (!frexpTest(FALSE)) {
  	warning("frexp: failed self-test")
  }
  if (!perturbTest()) {
  	warning("perturb: failed self-test")
  } 
  if (!secholTest()) {
  	warning("sechol: failed self-test")
  }
  return(TRUE)
}
