"aicc" <- function(model)
{
# Calculates AICc for a nls model fit (e.g. pw.sharp)
# Reference: Anderson, D.R., K.P. Burnham and W.L. Thompson. 2002. Null
# hypothesis testing: problems, prevalence, and an alternative. Journal of Wildlife
# Management 64:912-923.
#
# model = nls model fit
#
	sumfit <- summary(model)
	N<-sum(sumfit$df) 
	aicc <- N * log(sumfit$df[2]*sumfit$sigma^2/N) + 2 * sumfit$df[1] + 2*
sumfit$df[1]*( sumfit$df[1]+1)/(N- sumfit$df[1]-1)
	 return(aicc)
}

"aicc.fixedgamma" <- function(model)
{
# Calculates AICc for a nls model fit (e.g. pw.sharp) where gamma was fixed.
# Thus need to adjust for one extra estimted parameter (change the error and
# parameter df)
# Reference: Anderson, D.R., K.P. Burnham and W.L. Thompson. 2002. Null
# hypothesis testing: problems, prevalence, and an alternative. Journal of Wildlife
# Management 64:912-923.
#
# model = nls model fit
#
	sumfit <- summary(model)
	N <- sum(sumfit$df)
       K<-sumfit$df[1]+1
	aicc <- N * log((sumfit$df[2] * sumfit$sigma^2)/N) + 2 * K + (2 * K * (K +
1))/(N - K - 1)
	return(aicc)
}
"basic.q" <- function(z, g)
{
# From Chiu 2002 Thesis
# one-parameter basic bent cable function
	if(g > 0) ((z + g)^2/(4 * g)) * (z >=  - g & z <= g) + z * (z > g) else z * (z > 0)
}

"break.boot.ci" <- function(B, x, y, fnls, p, wts, ...)
{
# Bootstrap the regression Data using the empirical distribution of the residuals.
# This function assumes that the variance of the errors is constant.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# B = number of "good" bootstrap samples required
# fnls = appropriate piecewise model function (e.g. pw.sharp)
# p = vector of starting value(s) required for the fnls function
# wts = vector of weights to use in fitting
#
# rs is the random seed.  This can be used to replicate results.
# The original fit is in the first column of bootb
# Each row of failboot gives a bootstrap sample of y which did not
#    converge to an answer when fitted with fnls.  These models can
#    be examined to diagnose problems in the bootstrapping.
#
	rs <- .Random.seed
	datadf <- data.frame(x = x, y = y, weights = wts)
	data.nls <- fnls(datadf$x, datadf$y, p, datadf$weights, ...)	#original fit
	e <- resid(data.nls)
	bootb <- data.frame(coef(data.nls))
	failboot <- NULL
	n <- length(datadf$y)
	yfit <- fitted(data.nls)
	d <- datadf
	while(ncol(bootb) < (B + 1)) {
		j <- sample(n, replace = T)	#resample residuals
		d$y <- (yfit + e[j])/sqrt(d$weights)	#add to fitted y to give a bootstrapped y
		boot.nls <- fnls(d$x, d$y, p, d$weights, ...)
		if(length(boot.nls) != 1)
			bootb <- cbind(bootb, coef(boot.nls))	
	#if doesn't fail to converge, save coefficients
		if(length(boot.nls) == 1)
			failboot <- rbind(failboot, d$y)	#if fails to converge, save bootstrapped y
	}
	return(rs, bootb, failboot)
}

"gen.lin.hyp" <- function(contrast, fit)
{
# Runs a t-test for the hypothesis that a general linear combination
#    of model parameters equals zero.
# contrast = vector representing the linear combination of parameters
#    e.g. c(0,0,1,1) represents the combination 0*b0 + 0*b1 + 1*b2 + 1*b3 = 0
# fit = model fitted from a nls model (e.g. pw.sharp)
#
	summ.fit <- summary(fit)
	cov.mat <- summ.fit$cov.unscaled * summ.fit$sigma^2
	se.comb <- sqrt(t(contrast) %*% cov.mat %*% contrast)
	tc <- t(contrast) %*% fit$param/se.comb
	pvalue <- 2 * (1 - pt(abs(tc), summ.fit$df[2]))
	return(tc, pvalue)
}

"invert.ftest.bentcab" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a bent cable
#    piecewise model.
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
#   as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
# Gamma is unconstrained
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# p0 = vecot of initial estimate of breakpoint and gamma; e.g. p0=c(br0,g0)
# weights = vector of weights to use in fitting
# startstep = number of se(br) from brmax in which to search for the root above
#   The nonlinear model is difficult to fit for fixed values of br far away from
#   the max.
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	INF <- 1000000000
	datadf <- data.frame(dy, dx, weights)
	nlsmax <- nls(sqrt(datadf$weights) * datadf$dy ~ sqrt(datadf$weights) * cbind(1, datadf$dx, basic.q(
		datadf$dx - tau, gamm)), start = list(tau = p0[1], gamm = p0[2]), algorithm = "plinear")
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["tau"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["tau", "tau"] * ds2)	
	# start at brmax-startstep*se(br) and search towards brmax for root
	n <- length(dy)
	const <-  - error.df - qf(alpha, 1, error.df)
	dxsort <- sort(unique(dx))
	br <- brmax - startstep * sebr
	sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx, basic.q(datadf$dx - 
		br, gamm)), start = list(gamm = p0[2]), algorithm = "plinear"))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brlower <- br
	else {
		br <- br + xeps
		while((sign(obj1) == sign(obj2)) & (br < brmax)) {
			obj1 <- obj2
			sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx, 
				basic.q(datadf$dx - br, gamm)), start = list(gamm = p0[2]), algorithm = 
				"plinear"))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br + xeps
		}
		if(br >= brmax)
			brlower <-  - INF
		else brlower <- br - 2 * xeps
	}
	br <- brmax + startstep * sebr
	sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx, basic.q(datadf$dx - 
		br, gamm)), start = list(gamm = p0[2]), algorithm = "plinear"))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brupper <- br
	else {
		br <- br - xeps
		while((sign(obj1) == sign(obj2)) & (br > brmax)) {
			obj1 <- obj2
			sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx, 
				basic.q(datadf$dx - br, gamm)), start = list(gamm = p0[2]), algorithm = 
				"plinear"))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br - xeps
		}
		if(br <= brmax)
			brupper <- INF
		else brupper <- br + 2 * xeps
	}
	return(brmax, sebr, brlower, brupper)
}

"invert.ftest.bentcab.x" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95)
{
#
# Special modification for the case where gamma is fixed. 
#
# Inverted F-test confidence interval for the breakpoint in a bent cable
#    piecewise model.
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
#   as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# p0 = vecot of initial estimate of breakpoint and gamma; e.g. p0=c(br0,g0)
# weights = vector of weights to use in fitting
# startstep = number of se(br) from brmax in which to search for the root above
#   The nonlinear model is difficult to fit for fixed values of br far away from
#   the max.
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	INF <- 1000000000
	datadf <- data.frame(dy, dx, weights, g = p0[2])
	g <- p0[2]
	nlsmax <- nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx, basic.q(dx - tau, g[1])), data = 
		datadf, start = list(tau = p[1]), algorithm = "plinear")
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["tau"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["tau", "tau"] * ds2)	
	# start at brmax-startstep*se(br) and search towards brmax for root
	n <- length(dy)
	const <-  - error.df - qf(alpha, 1, error.df)
	dxsort <- sort(unique(dx))
	br <- brmax - startstep * sebr
	sum.fit <- summary(lm(datadf$dy ~ cbind(dx, basic.q(datadf$dx - br, g)), weights = datadf$weights))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brlower <- br
	else {
		br <- br + xeps
		while((sign(obj1) == sign(obj2)) & (br < brmax)) {
			obj1 <- obj2
			sum.fit <- summary(lm(datadf$dy ~ cbind(dx, basic.q(datadf$dx - br, g)), weights = 
				datadf$weights))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br + xeps
		}
		if(br >= brmax)
			brlower <-  - INF
		else brlower <- br - 2 * xeps
	}
	br <- brmax + startstep * sebr
	sum.fit <- summary(lm(datadf$dy ~ cbind(dx, basic.q(datadf$dx - br, g)), weights = datadf$weights))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brupper <- br
	else {
		br <- br - xeps
		while((sign(obj1) == sign(obj2)) & (br > brmax)) {
			obj1 <- obj2
			sum.fit <- summary(lm(datadf$dy ~ cbind(dx, basic.q(datadf$dx - br, g)), weights = 
				datadf$weights))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br - xeps
		}
		if(br <= brmax)
			brupper <- INF
		else brupper <- br + 2 * xeps
	}
	return(brmax, sebr, brlower, brupper)
}

"invert.ftest.benthyp" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95)
{
#
# find L.S. estimates of b0, b1, b2 and br in the bent hyperbola piecewise model
# dy = b0 + b1*(dx-br) + b2*sqrt((dx-br)^2+g^2/4)
# p0 = vector of initial estimates of br and g;  e.g. p0=c(br0,g0)
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
#   as alpha% confidence interval for br (i.e. alpha=.95 for 95% c.i.)
# Gamma is unconstrained
# S(beta(br), br) is the constrained SSerror when br has value br
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# startstep = number of se(br) from brmax in which to search for the root above
#   The nonlinear model is difficult to fit for fixed values of br far away from
#   the max.
# sample call:  
#
	INF <- 1000000000
	datadf <- data.frame(dy, dx, weights)
	nlsmax <- pw.benthyp(datadf$dx, datadf$dy, p0, datadf$weights)
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)	
	# start at brmax-startstep*se(br) and search towards brmax for root
	n <- length(dy)
	const <-  - error.df - qf(alpha, 1, error.df)
	dxsort <- sort(unique(dx))
	br <- brmax - startstep * sebr
	sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br, sqrt((distance -
		br)^2 + g^2/4)), start = list(g = p0[2]), algorithm = "plinear"))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brlower <- br
	else {
		br <- br + xeps
		while((sign(obj1) == sign(obj2)) & (br < brmax)) {
			obj1 <- obj2
			sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br,
				sqrt((distance - br)^2 + g^2/4)), start = list(g = p0[2]), algorithm = 
				"plinear"))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br + xeps
		}
		if(br >= brmax)
			brlower <-  - INF
		else brlower <- br - 2 * xeps
	}
	br <- brmax + startstep * sebr
	sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br, sqrt((distance -
		br)^2 + g^2/4)), start = list(g = p0[2]), algorithm = "plinear"))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brupper <- br
	else {
		br <- br - xeps
		while((sign(obj1) == sign(obj2)) & (br > brmax)) {
			obj1 <- obj2
			sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br,
				sqrt((distance - br)^2 + g^2/4)), start = list(g = p0[2]), algorithm = 
				"plinear"))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br - xeps
		}
		if(br <= brmax)
			brupper <- INF
		else brupper <- br + 2 * xeps
	}
	return(brmax, sebr, brlower, brupper)
}

"invert.ftest.sharp" <- function(dx, dy, br0, weights, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a sharp piecewise model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
#   as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when br has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	INF <- 1000000000
	datadf <- data.frame(dy, dx, weights)
	nlsmax <- pw.sharp(datadf$dx, datadf$dy, br0, datadf$weights)
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)	
	# need to start at the 2nd smallest observed dx and search for root
	n <- length(dy)
	const <-  - error.df - qf(alpha, 1, error.df)
	dxsort <- sort(unique(dx))
	br <- dxsort[2]
	sum.fit <- summary(lm(datadf$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights = weights))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brlower <- br
	else {
		br <- br + xeps
		while((sign(obj1) == sign(obj2)) & (br < brmax)) {
			obj1 <- obj2
			sum.fit <- summary(lm(datadf$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights
				 = weights))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br + xeps
		}
		if(br >= brmax)
			brlower <-  - INF
		else brlower <- br - 2 * xeps
	}
	br <- rev(dxsort)[2]
	sum.fit <- summary(lm(datadf$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights = weights))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brupper <- br
	else {
		br <- br - xeps
		while((sign(obj1) == sign(obj2)) & (br > brmax)) {
			obj1 <- obj2
			sum.fit <- summary(lm(datadf$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights
				 = weights))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br - xeps
		}
		if(br <= brmax)
			brupper <- INF
		else brupper <- br + 2 * xeps
	}
	return(brmax, sebr, brlower, brupper)
}

"invert.ftest.tanh" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a hyperbolic tangent
#    piecewise model.
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
#   as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
# Gamma is unconstrained
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# p0 = vecot of initial estimate of breakpoint and gamma; e.g. p0=c(br0,g0)
# weights = vector of weights to use in fitting
# startstep = number of se(br) from brmax in which to search for the root above
#   The nonlinear model is difficult to fit for fixed values of br far away from
#   the max.
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	INF <- 1000000000
	datadf <- data.frame(dy, dx, weights)
	nlsmax <- pw.tanh(datadf$dx, datadf$dy, p0, datadf$weights)
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)	
	# start at brmax-startstep*se(br) and search towards brmax for root
	n <- length(dy)
	const <-  - error.df - qf(alpha, 1, error.df)
	dxsort <- sort(unique(dx))
	br <- brmax - startstep * sebr
	sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br, (dx - br) * 
		tanh((dx - br)/g)), start = list(g = p0[2]), algorithm = "plinear"))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brlower <- br
	else {
		br <- br + xeps
		while((sign(obj1) == sign(obj2)) & (br < brmax)) {
			obj1 <- obj2
			sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br,
				(dx - br) * tanh((dx - br)/g)), start = list(g = p0[2]), algorithm = 
				"plinear"))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br + xeps
		}
		if(br >= brmax)
			brlower <-  - INF
		else brlower <- br - 2 * xeps
	}
	br <- brmax + startstep * sebr
	sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br, (dx - br) * 
		tanh((dx - br)/g)), start = list(g = p0[2]), algorithm = "plinear"))
	obj1 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
	obj2 <- obj1
	if(obj1 == 0)
		brupper <- br
	else {
		br <- br - xeps
		while((sign(obj1) == sign(obj2)) & (br > brmax)) {
			obj1 <- obj2
			sum.fit <- summary(nls(sqrt(weights) * datadf$dy ~ sqrt(weights) * cbind(1, dx - br,
				(dx - br) * tanh((dx - br)/g)), start = list(g = p0[2]), algorithm = 
				"plinear"))
			obj2 <- (sum.fit$sigma^2 * sum.fit$df[2])/ds2 + const
			br <- br - xeps
		}
		if(br <= brmax)
			brupper <- INF
		else brupper <- br + 2 * xeps
	}
	return(brmax, sebr, brlower, brupper)
}

"lack.of.fit" <- function(x, y, fit, weights)
{
# tests of lack-of-fit for nonlinear models versus saturated model
#
# x = vector of x data to fit
# y = vector of y data to fit
# fit = model fit to be tested fot lack of fit (e.g. resulting from pw.sharp)
# weights = vector of weights to use in fitting
#
	factor.fit <- lm(y ~ factor(x), weights = weights)
	dfpe <- factor.fit$df.residual
	sspe <- (summary(factor.fit))$sigma^2 * dfpe
	sse <- (summary(fit))$sigma^2 * (summary(fit))$df[2]
	dflof <- (summary(fit))$df[2] - dfpe
	Fc <- ((sse - sspe) * dfpe)/dflof/sspe
	pvalue <- 1 - pf(Fc, dflof, dfpe)
	return(Fc, pvalue)
}

"llsurface.bentcab" <- function(dx, dy, br, g, startvals, weights)
{
# Calculates a linear fit for each grid point (grid of possible br values).
# Calculates a log-likelihood surface across the grid.
# For the bent cable piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
#    persp(out$br,out$gamma,out$loglik,xlab="Value of the breakpoint, br", 
#       ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting values for breakpoint and gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting
#
	y <- dy
	x <- dx
	datadf <- data.frame(x = dx, y = dy, weights = weights)
	p0 <- startvals
	nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, distance, basic.q(distance - br, g)), 
		start = list(br = p0[1], g = p0[2]), data = datadf, algorithm = "plinear")
	nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]	#sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
#n <- length(y)	#const <- 4. - n - qf(alpha, 1., n - 4.)
	out <- NA
	out$loglik <- matrix(nrow = length(br), ncol = length(g))	
	# walk through all values of br and gamma, calculating the deviance
	for(i in 1:length(br)) {
		for(j in 1:length(g)) {
			out$loglik[i, j] <- summary(lm(datadf$y ~ cbind(distance, basic.q(distance - br[i], 
				g[j])), weights = weights))$sigma^2
		}
	}
	out$loglik <- (out$loglik * (error.df + 2) - ds2 * error.df)/ds2/2
	out$br <- br
	out$gamma <- g
	return(out)
}

"llsurface.bentcab.x" <- function(dx, dy, br, g, startvals, weights)
{
#
# Special modification for the case where gamma is fixed.
#
# Calculates a linear fit for each grid point (grid of possible br values).
# Calculates a log-likelihood surface across the grid.
# For the bent cable piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
#    persp(out$br,out$gamma,out$loglik,xlab="Value of the breakpoint, br", 
#       ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting value for breakpoint, fixed value for gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting
#
	y <- dy
	x <- dx
	p0 <- startvals
	g0 <- p0[2]	#datadf <- data.frame(x = dx, y = dy, weights = weights, p0 = startvals)
	datadf <- data.frame(x = dx, y = dy, weights = weights, g0)	
	#nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1., x, basic.q(x - br, p0[2.])), start #= list(br = p0[1.]), data = datadf, algorithm = "plinear")
	nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, basic.q(x - br, g0[1])), start = list(
		br = p0[1]), data = datadf, algorithm = "plinear")
	nlsum <- summary(nlsmax)
	ds2 <- nlsum$sigma^2	#sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
#n <- length(y) #const <- 4. - n - qf(alpha, 1., n - 4.)
	error.df <- nlsum$df[2]
	out <- NA	# walk through all values of br and gamma, calculating the deviance
	out$loglik <- matrix(nrow = length(br), ncol = length(g))
	for(i in 1:length(br)) {
		for(j in 1:length(g)) {
			out$loglik[i, j] <- summary(lm(datadf$y ~ cbind(datadf$x, basic.q(datadf$x - br[i], 
				g[j])), weights = weights))$sigma^2
		}
	}
	out$loglik <- (out$loglik * (error.df + 1) - ds2 * error.df)/ds2/2
	out$br <- br
	out$gamma <- g
	return(out)
}

"llsurface.benthyp" <- function(dx, dy, br, g, startvals, weights)
{
# Calculates a linear fit for each grid point (grid of possible br values).
# Plots a log-likelihood surface across the grid.
# For the bent hyperbola piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
#    persp(out$br,out$gamma,out$loglik,xlab="Value of the breakpoint, br", 
#       ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting values for breakpoint and gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting
#
	y <- dy
	x <- dx
	datadf <- data.frame(x = dx, y = dy, weights = weights)
	p0 <- startvals
	nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), sqrt((distance - br)^2 + 
		g^2/4)), start = list(br = p0[1], g = p0[2]), data = datadf, algorithm = "plinear")
	nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]	#sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
#n <- length(y)	#const <- 4. - n - qf(alpha, 1., n - 4.)
	out <- NA
	out$loglik <- matrix(nrow = length(br), ncol = length(g))	
	# walk through all values of br and gamma, calculating the deviance
	for(i in 1:length(br)) {
		for(j in 1:length(g)) {
			out$loglik[i, j] <- summary(lm(datadf$y ~ cbind((distance - br[i]), sqrt((distance - 
				br[i])^2 + g^2/4)), weights = weights))$sigma^2
		}
	}
	out$loglik <- (out$loglik * (error.df + 2) - ds2 * error.df)/ds2/2
	out$br <- br
	out$gamma <- g
	return(out)
}

"llsurface.sharp" <- function(dx, dy, br, br.start, weights)
{
# Calculates a linear fit for each grid point (grid of possible br values).
# Calculates a log-likelihood surface across the grid.
# For the sharp piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
#    persp(out$br,out$gamma,out$loglik,xlab="Value of the breakpoint, br", 
#       ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# br.start = initial starting values for breakpoint; e.g. br.start=br0
# weights = vector of weights to use in fitting
#
	y <- dy
	x <- dx
	datadf <- data.frame(x = dx, y = dy, weights = weights)
	br0 <- br.start
	nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, ifelse(x > br, x - br, 0)), start = 
		list(br = br0), control = list(maxiter = 15, tolerance = 1e-005, minscale = 1e-005), data
		 = datadf, algorithm = "plinear")
	nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]	#sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
#n <- length(y) #const <- 4. - n - qf(alpha, 1., n - 4.)
	out <- data.frame(br = br, loglik = rep(NA, length(br)))	
	# walk through all values of br, calculating the deviance
	out$loglik <- sapply(br, loopfn.sharp, datadf, weights)
	out$loglik <- (out$loglik * (error.df + 1) - ds2 * error.df)/ds2/2
	print(cbind(br.fitted = brmax, sebr.fitted = sebr, sigmasq.fitted = ds2))
	return(out)
}

"llsurface.tanh" <- function(dx, dy, br, g, startvals, weights)
{
# Calculates a linear fit for each grid point (grid of possible br, gamma)
# Plots a scaled deviance surface across the grid, 
# for the hyperbolic tangent piecewise regression model.
# This surface can be used to compute approx. conf regions by comparing with
# quantiles of the F(2, nlsum$df[2]), since the scaling incorporates the 2 df
# for the 2 parameters br and gamma.
#
# Can use to create a perspective log-likelihood plot:
#    persp(out$br,out$gamma,out$loglik,xlab="Value of the breakpoint, br", 
#       ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting values for breakpoint and gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting; proportional to 1/Var(dy)
#
	y <- dy
	x <- dx
	datadf <- data.frame(x = dx, y = dy, weights = weights)
	p0 <- startvals
	nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), (distance - br) * tanh((
		distance - br)/g)), start = list(br = p0[1], g = p0[2]), data = datadf, algorithm = 
		"plinear")
	nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]	#sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
#n <- length(y)	#const <- 4. - n - qf(alpha, 1., n - 4.)
	out <- NA
	out$loglik <- matrix(nrow = length(br), ncol = length(g))	
	# walk through all values of br and gamma, calculating the deviance
	for(i in 1:length(br)) {
		for(j in 1:length(g)) {
			out$loglik[i, j] <- summary(lm(datadf$y ~ cbind((distance - br[i]), (distance - br[
				i]) * tanh((distance - br[i])/g[j])), weights = weights))$sigma^2
		}
	}
	out$loglik <- (out$loglik * (error.df + 2) - ds2 * error.df)/ds2/2
	out$br <- br
	out$gamma <- g
	return(out)
}

"loopfn.sharp" <- function(br, data, weights)
{
# internal function for llsurface.sharp
#
	datadf <- data
	br0 <- br
	loglik <- summary(lm(datadf$y ~ cbind(x, ifelse(x > br0, x - br0, 0)), weights = weights))$sigma^2
	return(loglik)
}

"mynls" <- function(formula, data = sys.parent(), start = if(!is.null(attr(data, "parameters"))) attr(data, 
		"parameters") else nls.initial(), control, algorithm = "default", trace = F)
{
	convert.twiddle <- function(formula)
	{
		if(length(formula) < 3)
			return(formula)
		form <- call("~", call("-", formula[[2]], formula[[3]]))
		attr(form, "class") <- "formula"
		form
	}
	if(is.numeric(data))
		data <- sys.frame(data)
	cl <- class(data)
	if(inherits(data, "pframe")) {
		class(data) <- NULL
		pp <- parameters(data)
		if(length(pp)) {
			np <- names(pp)
			if(any(match(np, names(data), 0)))
				stop("can't have variables, parameters with\nsame na\nme")
			data[np] <- pp
		}
	}
	else if(inherits(data, "data.frame"))
		cl <- c("pframe", cl)
	class(data) <- NULL	## First, figure out the par. names, make start a list
	switch(mode(start),
		numeric = {
			.parameters <- names(start)
			start <- as.list(start)
			names(start) <- .parameters
		}
		,
		list = {
			.parameters <- names(start)
		}
		,
		NULL = .parameters <- parameter.names(formula, data),
		stop("\"start\" should be numeric or list"))
	if(!length(.parameters))
		stop("names for parameters needed, from formula or from start")
	pn <- .parameters
	.expr <- formula	## select the algorithm and possibly massage the formula
	if(length(start)) data[.parameters] <- start	#used by setup_nonlin
	data$.parameters <- .parameters
	nl.frame <- new.frame(data, F)
	frame.attr("class", nl.frame) <- cl	# in case data is returned
	formula <- switch(algorithm,
		plinear = {
			if(length(formula) < 3)
				stop("formula for plinear algorithm be of the form\nresp ~ array")
			response <- eval(formula[[2]], nl.frame)
			design <- eval(formula[[3]], nl.frame)
			nnobs <- length(response)
			nlinear <- if(is.matrix(design)) dim(design)[2] else length(design)/nnobs
			form <- call("~", formula[[3]])
			attr(form, "class") <- "formula"
			form
		}
		,
		default = convert.twiddle(formula))
	dims <- .C("setup_nonlin",
		n = integer(3),
		list(formula),
		as.integer(nl.frame))$n
	npar <- dims[1]
	nderiv <- dims[2]
	nobs <- dims[3]
	resp <- eval(.expr[[2]], nl.frame)
	if(length(.expr) == 2) resp[] <- 0	## setup_nonlin will have set .parameters if missing
	if(is.null(start)) {
		start <- list()
		if(is.null(names(start)) && length(start) == length(pn))
			names(start) <- pn
		for(i in .parameters)
			start[[i]] <- get(i, frame = nl.frame)
	}
	asgn <- start
	si <- 1
	for(i in 1:length(asgn)) {
		ni <- length(asgn[[i]])
		asgn[[i]] <- seq(from = si, length = ni)
		si <- si + ni
	}
	start <- unlist(start)
	controlvals <- nls.control()
	if(!missing(control))
		controlvals[names(control)] <- control
	ret.data <- is.logical(ret.data <- controlvals$data) && ret.data
	max.iterations <- if(is.null(controlvals$maxiter)) 50 * npar else controlvals$maxiter
	settings <- c(0, max.iterations, controlvals$minscale, controlvals$tolerance, 0)
	if(algorithm == "plinear") {
		settings[1] <- 1
		dims <- c(dims, nlinear)
		start <- c(start, numeric(nlinear))
		pn <- c(pn, paste(".lin", 1:nlinear, sep = ""))
		dims[3] <- nnobs
		outmat <- array(c(response, numeric(nnobs * (npar + nlinear))), c(nnobs, npar + nlinear + 1
			))
	}
	else outmat <- array(0, c(nobs, npar + 1))
	storage.mode(outmat) <- "double"
	storage.mode(start) <- "double"
	nls.trace <- if(missing(trace)) controlvals$trace else trace
	std.trace <- FALSE
	if(is.logical(nls.trace)) {
		if(std.trace <- nls.trace)
			nls.trace <- "trace.nls"
		else nls.trace <- NULL
	}
	else std.trace <- is.character(nls.trace) && nls.trace == "trace.nls"
	if(std.trace) {
		assign("trace.mat", array(0, c(max.iterations, npar + 2), list(NULL, c("obj.", "conv.", 
			paste("par", 1:npar)))), frame = 1)
		assign("trace.expr", expression(trace.mat[last.iteration,  ] <- it.row), frame = 1)
	}
	z <- .C("do_nls",
		parameters = start,
		dims = as.integer(dims),
		control = as.double(settings),
		outmat = outmat,
		trace = list(nls.trace))
	if(z$control[5]) nls.out <- c("skip") else {
		nls.out <- list(parameters = z$parameters, formula = .expr, call = match.call(), residuals
			 = z$outmat[, 1])
		if(algorithm == "plinear") {
			class(nls.out) <- c("nls.pl", "nls")
			R <- qr(z$outmat[, -1])$qr[1:(npar + nlinear),  , drop = F]
		}
		else {
			class(nls.out) <- "nls"
			R <- qr(z$outmat[, -1])$qr[1:npar,  , drop = F]
		}
		R[lower.tri(R)] <- 0
		nls.out$R <- R
		nls.out$fitted.values <- resp - nls.out$residuals
		nls.out$assign <- asgn
		if(ret.data) {
			data <- sys.frame(nl.frame)
			data$.parameters <- NULL	
	#dbdetach does the right thing--only matters that nl.frame>1
			if(inherits(data, "pframe"))
				data <- dbdetach(data, nl.frame)
			nls.out$data <- data
		}
		if(std.trace && exists("last.iteration", frame = 1))
			nls.out$trace <- get("trace.mat", frame = 1)[1:get("last.iteration", frame = 1),  ]
	}
#       stop(switch(as.integer(z$control[5]), 
#               "step factor reduced below minimum", 
#               "maximum number of iterations exceeded", 
#               "singular gradient matrix", 
#               "singular design matrix", 
#               "singular gradient matrix"))
	nls.out
}

"plot.invert.ftest.bentcab" <- function(dx, dy, br0, g0, weights, limits, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a bent cable piecewise
#    model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
#    as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# limits = vector length 2 of start and end limits
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	datadf <- data.frame(dy, dx, weights, br = br0)
	p0 <- c(br0, g0)
	nlsmax <- pw.bentcab(datadf$dx, datadf$dy, p0, datadf$weights)
	print(nlsmax)
	nlsum <- summary(nlsmax)
	print(nlsum)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
	print(brmax)
	Fcrit <- qf(alpha, 1, error.df)	# undx <- sort(unique(dx))
# brgrid <- seq(undx[2], rev(undx)[2], xeps)
	brgrid <- seq(limits[1], limits[2], xeps)
	obj1 <- NA
	for(br in brgrid) {
		print(br)
		datadf$br <- br
		sum.fit <- summary(nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx, basic.q(dx - br, 
			g0)), start = list(g0 = p0[2]), data = datadf, algorithm = "plinear"))
		obj1 <- cbind(obj1, sum.fit$sigma)
	}
	Fstats <- (obj1[-1]^2 * sum.fit$df[2])/ds2 - error.df
	Fstats <- cbind(brgrid, Fstats)
	plot(Fstats[, 1], Fstats[, 2], type = "l")
	abline(h = Fcrit)
	return(brmax, sebr, Fcrit, Fstats)
}

"plot.invert.ftest.bentcab.x" <- function(dx, dy, br0, g0, weights, xeps = 2, alpha = 0.95)
{
# Special modification for the case where gamma is fixed.
#
# Inverted F-test confidence interval for the breakpoint in a bent cable piecewise
#    model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
#    as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	p0 <- c(br0, g0)
	datadf <- data.frame(dy, dx, weights, p0)
	nlsmax <- nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx, basic.q(dx - br, p0[2])), data = 
		datadf, start = list(br = p0[1]), algorithm = "plinear")
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
	Fcrit <- qf(alpha, 1, error.df)
	undx <- sort(unique(dx))
	brgrid <- seq(undx[2], rev(undx)[2], xeps)
	obj1 <- NA
	for(br in brgrid) {
		sum.fit <- summary(lm(datadf$dy ~ cbind(distance, basic.q(distance - br, p0[2])), weights
			 = weights))
		obj1 <- cbind(obj1, sum.fit$sigma)
	}
	Fstats <- (obj1[-1]^2 * sum.fit$df[2])/ds2 - error.df
	Fstats <- cbind(brgrid, Fstats)
	plot(Fstats[, 1], Fstats[, 2], type = "l")
	abline(h = Fcrit)
	return(brmax, sebr, Fcrit, Fstats)
}

"plot.invert.ftest.benthyp" <- function(dx, dy, br0, g0, weights, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a bent hyperbola piecewise
#    model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
#    as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	datadf <- data.frame(dy, dx, weights)
	p0 <- c(br0, g0)
	nlsmax <- pw.benthyp(datadf$dx, datadf$dy, p0, datadf$weights)
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
	Fcrit <- qf(alpha, 1, error.df)
	undx <- sort(unique(dx))
	brgrid <- seq(undx[2], rev(undx)[2], xeps)
	obj1 <- NA
	for(br in brgrid) {
		sum.fit <- summary(lm(datadf$dy ~ cbind(distance - br, sqrt((distance - br)^2 + (g0)^2/4)), 
			weights = weights))
		obj1 <- cbind(obj1, sum.fit$sigma)
	}
	Fstats <- (obj1[-1]^2 * sum.fit$df[2])/ds2 - error.df
	Fstats <- cbind(brgrid, Fstats)
	plot(Fstats[, 1], Fstats[, 2], type = "l")
	abline(h = Fcrit)
	return(brmax, sebr, Fcrit, Fstats)
}

"plot.invert.ftest.sharp" <- function(dx, dy, br0, weights, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a sharp piecewise
#    model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
#    as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but 
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	datadf <- data.frame(dy, dx, weights)
	nlsmax <- pw.sharp(datadf$dx, datadf$dy, br0, datadf$weights)
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
	Fcrit <- qf(alpha, 1, error.df)
	undx <- sort(unique(dx))
	brgrid <- seq(undx[2], rev(undx)[2], xeps)
	obj1 <- NA
	for(br in brgrid) {
		sum.fit <- summary(lm(datadf$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights = weights
			))
		obj1 <- cbind(obj1, sum.fit$sigma)
	}
	Fstats <- (obj1[-1]^2 * sum.fit$df[2])/ds2 - error.df
	Fstats <- cbind(brgrid, Fstats)
	plot(Fstats[, 1], Fstats[, 2], type = "l")
	abline(h = Fcrit)
	return(brmax, sebr, Fcrit, Fstats)
}

"plot.invert.ftest.tanh" <- function(dx, dy, br0, g0, weights, limits, xeps = 2, alpha = 0.95)
{
# Inverted F-test confidence interval for the breakpoint in a hyperbolic
#  tangent piecewise
#    model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2
#  <=F(alpha,1,error.df)
#    as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# limits = vector length 2 of start and end limits
# xeps = step size in grid search, smaller values are more precise but
#   computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#
	datadf <- data.frame(dy, dx, weights, br = br0)
	p0 <- c(br0, g0)
	nlsmax <- pw.tanh(datadf$dx, datadf$dy, p0, datadf$weights)
	nlsum <- summary(nlsmax)
	brmax <- (coef(nlsmax))["br"]
	ds2 <- nlsum$sigma^2
	error.df <- nlsum$df[2]
	sebr <- sqrt(nlsum$cov.unscaled["br", "br"] * ds2)
	print(brmax)
	Fcrit <- qf(alpha, 1, error.df)	#       undx <- sort(unique(dx))
#       brgrid <- seq(undx[2], rev(undx)[2], xeps)
	brgrid <- seq(limits[1], limits[2], xeps)
	obj1 <- NA
	for(br in brgrid) {
		print(br)
		datadf$br <- br
		sum.fit <- summary(nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx - br, (dx - br) * 
			tanh((dx - br)/g0)), start = list(g0 = p0[2]), data = datadf, algorithm = "plinear"
			))
		obj1 <- cbind(obj1, sum.fit$sigma)
	}
	Fstats <- (obj1[-1]^2 * sum.fit$df[2])/ds2 - error.df
	Fstats <- cbind(brgrid, Fstats)
	plot(Fstats[, 1], Fstats[, 2], type = "l")
	abline(h = Fcrit)
	return(brmax, sebr, Fcrit, Fstats)
}

"pw.bentcab" <- function(x, y, p, weights, maxit = 50)
{
# Fits bent cable piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting values (breakpoint, gamma)
# weights = vector of weights to use in fitting
#
	datadf <- data.frame(x = x, y = y, weights = weights)
	data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, basic.q(x - br, g)), data = 
		datadf, start = list(br = p[1], g = p[2]), algorithm = "plinear", control = list(maxiter = 
		maxit))
	return(data.nls)
}

"pw.bentcab.x" <- function(x, y, p, weights, maxit = 50)
{
# Fits bent cable piecewise regression with one breakpoint.
# Modification for the case when gamma needs to be fixed.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting value of breakpoint, fixed value for gamma
# weights = vector of weights to use in fitting
#
	datadf <- data.frame(x = x, y = y, weights = weights, g = p[2])
	data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, basic.q(x - br, g[1])), data = 
		datadf, start = list(br = p[1]), algorithm = "plinear", control = list(maxiter = maxit))
	return(data.nls)
}

"pw.benthyp" <- function(x, y, p, weights, maxit = 50)
{
# Fits bent hyperbola piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting values (breakpoint, gamma)
# weights = vector of weights to use in fitting
#
	datadf <- data.frame(x = x, y = y, weights = weights)
	data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), sqrt((distance - br
		)^2 + g^2/4)), data = datadf, start = list(br = p[1], g = p[2]), algorithm = "plinear", 
		control = list(maxiter = maxit))
	return(data.nls)
}

"pw.sharp" <- function(x, y, p, weights, minscale = 0.001, toler = 0.001, trace = F)
{
# Fits sharp piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = starting value of breakpoint
# weights = vector of weights to use in fitting
#
	datadf <- data.frame(x = x, y = y, weights = weights)
	data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, distance, ifelse(distance > br, 
		distance - br, 0)), data = datadf, start = list(br = p), algorithm = "plinear", control = c(
		minscale = minscale, tolerance = toler), trace = trace)
	return(data.nls)
}

"pw.tanh" <- function(x, y, p, weights, maxit = 50)
{
# Fits hyperbolic tangent piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting values (breakpoint, gamma)
# weights = vector of weights to use in fitting
#
	datadf <- data.frame(x = x, y = y, weights = weights)
	data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), (distance - br) * 
		tanh((distance - br)/g)), data = datadf, start = list(br = p[1], g = p[2]), algorithm = 
		"plinear", control = list(maxiter = maxit))
	return(data.nls)
}