#' Functions for positioning tick labels on axes
#'
#' \tabular{ll}{
#' Package: \tab labeling\cr
#' Type: \tab Package\cr
#' Version: \tab 0.2\cr
#' Date: \tab 2011-04-01\cr
#' License: \tab Unlimited\cr
#' LazyLoad: \tab yes\cr
#' }
#' 
#' Implements a number of axis labeling schemes, including those 
#' compared in An Extension of Wilkinson's Algorithm for Positioning Tick Labels on Axes
#' by Talbot, Lin, and Hanrahan, InfoVis 2010.
#'
#' @name labeling-package
#' @aliases labeling
#' @docType package
#' @title Axis labeling
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @references
#' Heckbert, P. S. (1990) Nice numbers for graph labels, Graphics Gems I, Academic Press Professional, Inc.
#' Wilkinson, L. (2005) The Grammar of Graphics, Springer-Verlag New York, Inc.
#' Talbot, J., Lin, S., Hanrahan, P. (2010) An Extension of Wilkinson's Algorithm for Positioning Tick Labels on Axes, InfoVis 2010.
#' @keywords dplot
#' @seealso \code{\link{extended}}, \code{\link{wilkinson}}, \code{\link{heckbert}}, \code{\link{rpretty}}, \code{\link{gnuplot}}, \code{\link{matplotlib}}, \code{\link{nelder}}, \code{\link{sparks}}, \code{\link{thayer}}, \code{\link{pretty}}
#' @examples
#' heckbert(8.1, 14.1, 4)	# 5 10 15
#' wilkinson(8.1, 14.1, 4)	# 8 9 10 11 12 13 14 15
#' extended(8.1, 14.1, 4)	# 8 10 12 14

#' # When plotting, extend the plot range to include the labeling
#' # Should probably have a helper function to make this easier
#' data(iris)
#' x <- iris$Sepal.Width
#' y <- iris$Sepal.Length
#' xl <- extended(min(x), max(x), 6)
#' yl <- extended(min(y), max(y), 6)
#' plot(x, y, 
#'     xlim=c(min(x,xl),max(x,xl)), 
#'     ylim=c(min(y,yl),max(y,yl)), 
#'     axes=FALSE, main="Extended labeling")
#' axis(1, at=xl)
#' axis(2, at=yl)
c()



#' Heckbert's labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @return vector of axis label locations
#' @references
#' Heckbert, P. S. (1990) Nice numbers for graph labels, Graphics Gems I, Academic Press Professional, Inc.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
heckbert <- function(dmin, dmax, m)
{
    range <- .heckbert.nicenum((dmax-dmin), FALSE)
    lstep <- .heckbert.nicenum(range/(m-1), TRUE)
    lmin <- floor(dmin/lstep)*lstep
    lmax <- ceiling(dmax/lstep)*lstep
    seq(lmin, lmax, by=lstep)
}
    
.heckbert.nicenum <- function(x, round)
{
	e <- floor(log10(x))
	f <- x / (10^e)
	if(round)
	{
		if(f < 1.5) nf <- 1
		else if(f < 3) nf <- 2
		else if(f < 7) nf <- 5
		else nf <- 10
	}
	else
	{
		if(f <= 1) nf <- 1
		else if(f <= 2) nf <- 2
		else if(f <= 5) nf <- 5
		else nf <- 10
	}
	nf * (10^e)
}



#' Wilkinson's labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @param Q set of nice numbers
#' @param mincoverage minimum ratio between the the data range and the labeling range, controlling the whitespace around the labeling (default = 0.8)
#' @param mrange range of \code{m}, the number of tick marks, that should be considered in the optimization search
#' @return vector of axis label locations
#' @note Ported from Wilkinson's Java implementation with some changes.
#'	Changes: 	1) m (the target number of ticks) is hard coded in Wilkinson's implementation as 5. 
#'						Here we allow it to vary as a parameter. Since m is fixed, 
#'						Wilkinson only searches over a fixed range 4-13 of possible resulting ticks.
#'						We broadened the search range to max(floor(m/2),2) to ceiling(6*m), 
#'						which is a larger range than Wilkinson considers for 5 and allows us to vary m,
#'						including using non-integer values of m.
#'				2) Wilkinson's implementation assumes that the scores are non-negative. But, his revised
#'						granularity function can be extremely negative. We tweaked the code to allow negative scores.
#'						We found that this produced better labelings.
#'				3) We added 10 to Q. This seemed to be necessary to get steps of size 1.
#'	It is possible for this algorithm to find no solution.
#'					In Wilkinson's implementation, instead of failing, he returns the non-nice labels spaced evenly from min to max.
#'					We want to detect this case, so we return NULL. If this happens, the search range, mrange, needs to be increased.
#' @references
#' Wilkinson, L. (2005) The Grammar of Graphics, Springer-Verlag New York, Inc.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
wilkinson <-function(dmin, dmax, m, Q = c(1,5,2,2.5,3,4,1.5,7,6,8,9), mincoverage = 0.8, mrange=max(floor(m/2),2):ceiling(6*m))
{
	best <- NULL
	for(k in mrange)
	{
		result <- .wilkinson.nice.scale(dmin, dmax, k, Q, mincoverage, mrange, m)
		if(!is.null(result) && (is.null(best) || result$score > best$score))
		{
			best <- result
		}
	}
	seq(best$lmin, best$lmax, by=best$lstep)
}

.wilkinson.nice.scale <- function(min, max, k, Q = c(1,5,2,2.5,3,4,1.5,7,6,8,9), mincoverage = 0.8, mrange=c(), m=k)
{
	Q <- c(10, Q)

	range <- max-min
	intervals <- k-1
	granularity <- 1 - abs(k-m)/m

	delta <- range / intervals
	base <- floor(log10(delta))
	dbase <- 10^base

	best <- NULL
	for(i in 1:length(Q))
	{
		tdelta <- Q[i] * dbase
		tmin <- floor(min/tdelta) * tdelta
		tmax <- tmin + intervals * tdelta

		if(tmin <= min && tmax >= max)
		{
			roundness <- 1 - ((i-1) - ifelse(tmin <= 0 && tmax >= 0, 1, 0)) / length(Q)
			coverage <- (max-min)/(tmax-tmin)
			if(coverage > mincoverage)
			{
				tnice <- granularity + roundness + coverage

				## Wilkinson's implementation contains code to favor certain ranges of labels
				## e.g. those balanced around or anchored at 0, etc.
				## We did not evaluate this type of optimization in the paper, so did not include it.
				## Obviously this optimization component could also be added to our function.
				#if(tmin == -tmax || tmin == 0 || tmax == 1 || tmax == 100)
				#	tnice <- tnice + 1
				#if(tmin == 0 && tmax == 1 || tmin == 0 && tmax == 100)
				#	tnice <- tnice + 1

				if(is.null(best) || tnice > best$score)
				{
					best <- list(lmin=tmin,
							 lmax=tmax,
							 lstep=tdelta,
					      	 score=tnice
						)
				}
			}
		}
	}
	best
}




## The Extended-Wilkinson algorithm described in the paper.

## Our scoring functions, including the approximations for limiting the search
.simplicity <- function(q, Q, j, lmin, lmax, lstep)
{
	eps <- .Machine$double.eps * 100

	n <- length(Q)
	i <- match(q, Q)[1]
	v <- ifelse( (lmin %% lstep < eps || lstep - (lmin %% lstep) < eps) && lmin <= 0 && lmax >=0, 1, 0)

	1 - (i-1)/(n-1) - j + v
}

.simplicity.max <- function(q, Q, j)
{
	n <- length(Q)
	i <- match(q, Q)[1]
	v <- 1

	1 - (i-1)/(n-1) - j + v
}

.coverage <- function(dmin, dmax, lmin, lmax)
{
	range <- dmax-dmin
	1 - 0.5 * ((dmax-lmax)^2+(dmin-lmin)^2) / ((0.1*range)^2)
}

.coverage.max <- function(dmin, dmax, span)
{
	range <- dmax-dmin
	if(span > range)
	{
		half <- (span-range)/2
		1 - 0.5 * (half^2 + half^2) / ((0.1 * range)^2)
	}
	else
	{
		1
	}
}

.density <- function(k, m, dmin, dmax, lmin, lmax)
{
	r <- (k-1) / (lmax-lmin)
	rt <- (m-1) / (max(lmax,dmax)-min(dmin,lmin))
	2 - max( r/rt, rt/r )
}

.density.max <- function(k, m)
{
	if(k >= m)
		2 - (k-1)/(m-1)
	else
		1
}

.legibility <- function(lmin, lmax, lstep)
{
	1			## did all the legibility tests in C#, not in R.
}


#' An Extension of Wilkinson's Algorithm for Position Tick Labels on Axes
#'
#' \code{extended} is an enhanced version of Wilkinson's optimization-based axis labeling approach. It is described in detail in our paper. See the references.
#' 
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @param Q set of nice numbers
#' @param only.loose if true, the extreme labels will be outside the data range
#' @param w weights applied to the four optimization components (simplicity, coverage, density, and legibility)
#' @return vector of axis label locations
#' @references
#' Talbot, J., Lin, S., Hanrahan, P. (2010) An Extension of Wilkinson's Algorithm for Positioning Tick Labels on Axes, InfoVis 2010.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
extended <- function(dmin, dmax, m, Q=c(1,5,2,2.5,4,3), only.loose=FALSE, w=c(0.25,0.2,0.5,0.05))
{
	eps <- .Machine$double.eps * 100
	
	if(dmin > dmax) {
		temp <- dmin
		dmin <- dmax
		dmax <- temp
	}

	if(dmax - dmin < eps) {
		#if the range is near the floating point limit,
		#let seq generate some equally spaced steps.
		return(seq(from=dmin, to=dmax, length.out=m))
	}
	
	if((dmax - dmin) > sqrt(.Machine$double.xmax)) {
    #if the range is too large
    #let seq generate some equally spaced steps.
    return(seq(from=dmin, to=dmax, length.out=m))
  }

	n <- length(Q)

	best <- list()
	best$score <- -2
	
	j <- 1
	while(j < Inf)
	{
		for(q in Q)
		{
			sm <- .simplicity.max(q, Q, j)

			if((w[1]*sm+w[2]+w[3]+w[4]) < best$score)
			{
				j <- Inf
				break
			}
		
			k <- 2
			while(k < Inf)													# loop over tick counts
			{		
				dm <- .density.max(k, m)

				if((w[1]*sm+w[2]+w[3]*dm+w[4]) < best$score)
					break
			
				delta <- (dmax-dmin)/(k+1)/j/q
				z <- ceiling(log(delta, base=10))

				while(z < Inf)
				{			
					step <- j*q*10^z

					cm <- .coverage.max(dmin, dmax, step*(k-1))

					if((w[1]*sm+w[2]*cm+w[3]*dm+w[4]) < best$score)
						break
					
					min_start <- floor(dmax/(step))*j - (k - 1)*j
					max_start <- ceiling(dmin/(step))*j

					if(min_start > max_start)
					{
						z <- z+1
						next
					}

					for(start in min_start:max_start)
					{
						lmin <- start * (step/j)
						lmax <- lmin + step*(k-1)
						lstep <- step

						s <- .simplicity(q, Q, j, lmin, lmax, lstep)
						c <- .coverage(dmin, dmax, lmin, lmax)						
						g <- .density(k, m, dmin, dmax, lmin, lmax)
						l <- .legibility(lmin, lmax, lstep)						

						score <- w[1]*s + w[2]*c + w[3]*g + w[4]*l

						if(score > best$score && (!only.loose || (lmin <= dmin && lmax >= dmax)))
						{
							best <- list(lmin=lmin,
								 lmax=lmax,
								 lstep=lstep,
						         score=score)
						}
					}
					z <- z+1
				}				
				k <- k+1
			}
		}
		j <- j + 1		
	}

	seq(from=best$lmin, to=best$lmax, by=best$lstep)
}



## Quantitative evaluation plots (Figures 2 and 3 in the paper)


#' Generate figures from An Extension of Wilkinson's Algorithm for Position Tick Labels on Axes
#'
#' Generates Figures 2 and 3 from our paper.
#' 
#' @param samples number of samples to use (in the paper we used 10000, but that takes awhile to run).
#' @return produces plots as a side effect
#' @references
#' Talbot, J., Lin, S., Hanrahan, P. (2010) An Extension of Wilkinson's Algorithm for Positioning Tick Labels on Axes, InfoVis 2010.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
extended.figures <- function(samples = 100)
{
	oldpar <- par()
	par(ask=TRUE)
	
	a <- runif(samples, -100, 400)
	b <- runif(samples, -100, 400)
	low <- pmin(a,b)
	high <- pmax(a,b)
	ticks <- runif(samples, 2, 10)

	generate.labelings <- function(labeler, dmin, dmax, ticks, ...)
	{
		mapply(labeler, dmin, dmax, ticks, SIMPLIFY=FALSE, MoreArgs=list(...))
	}
	
	h1 <- generate.labelings(heckbert, low, high, ticks)
	w1 <- generate.labelings(wilkinson, low, high, ticks, mincoverage=0.8)
	f1 <- generate.labelings(extended, low, high, ticks, only.loose=TRUE)
	e1 <- generate.labelings(extended, low, high, ticks)
	
	figure2 <- function(r, names)
	{
		for(i in 1:length(r))
		{
			d <- r[[i]]
			
			#plot coverage
			cover <- sapply(d, function(x) {max(x)-min(x)})/(high-low)
			hist(cover, breaks=seq(from=-0.01,to=1000,by=0.02), xlab="", ylab=names[i], main=ifelse(i==1, "Density", ""), col="darkgray", lab=c(3,3,3), xlim=c(0.5,3.5), ylim=c(0,0.12*samples), axes=FALSE, border=FALSE)
			#hist(cover)
			axis(side=1, at=c(0,1,2,3,4), xlab="hello", line=-0.1, lwd=0.5)
			
			# plot density
			dens <- sapply(d, length) / ticks
			hist(dens, breaks=seq(from=-0.01,to=10,by=0.02), xlab="", ylab=names[i], main=ifelse(i==1, "Density", ""), col="darkgray", lab=c(3,3,3), xlim=c(0.5,3.5), ylim=c(0,0.06*samples), axes=FALSE, border=FALSE)
			axis(side=1, at=c(0,1,2,3,4), xlab="hello", line=-0.1, lwd=0.5)
		}
	}

	par(mfrow=c(4, 2), mar=c(0.5,1.85,1,0), oma=c(1,0,1,0), mgp=c(0,0.5,-0.3), font.main=1, font.lab=1, cex.lab=1, cex.main=1, tcl=-0.2)
	figure2(list(h1,w1, f1, e1), names=c("Heckbert", "Wilkinson", "Extended\n(loose)", "Extended\n(flexible)"))

	figure3 <- function(r, names)
	{
		for(i in 1:length(r))
		{
			d <- r[[i]]
			steps <- sapply(d, function(x) round(median(diff(x)), 2))
			steps <- steps / (10^floor(log10(steps)))
			tab <- table(steps)
			barplot(rev(tab), xlim=c(0,0.4*samples), horiz=TRUE, xlab=ifelse(i==1,"Frequency",""), xaxt='n', yaxt='s', las=1, main=names[i], border=NA, col="gray")
		}
	}
	
	par(mfrow=c(1,4), mar=c(0.5, 0.75, 2, 0.5), oma=c(0,2,1,1), mgp=c(0,0.75,-0.3), cex.lab=1, cex.main=1)
	figure3(list(h1,w1, f1, e1), names=c("Heckbert", "Wilkinson", "Extended\n(loose)", "Extended\n(flexible)"))
	par(oldpar)
}



#' Nelder's labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @param Q set of nice numbers
#' @return vector of axis label locations
#' @references
#' Nelder, J. A. (1976) AS 96. A Simple Algorithm for Scaling Graphs, Journal of the Royal Statistical Society. Series C., pp. 94-96.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
nelder <- function(dmin, dmax, m, Q = c(1,1.2,1.6,2,2.5,3,4,5,6,8,10))
{
	ntick <- floor(m)
	tol <- 5e-6
	bias <- 1e-4

	intervals <- m-1
	x <- abs(dmax)
	if(x == 0) x <- 1
	if(!((dmax-dmin)/x > tol))
	{
		## special case handling for very small ranges. Not implemented yet.
	}

	step <- (dmax-dmin)/intervals
	s <- step

	while(s <= 1)
		s <- s*10
	while(s > 10)
		s <- s/10

	x <- s-bias
	unit <- 1
	for(i in 1:length(Q))
	{
		if(x < Q[i])
		{
			unit <- i
			break
		}
	}
	step <- step * Q[unit] / s
	range <- step*intervals

	x <- 0.5 * (1+ (dmin+dmax-range) / step)
	j <- floor(x-bias)
	valmin <- step * j

	if(dmin > 0 && range >= dmax)
		valmin <- 0
	valmax <- valmin + range

	if(!(dmax > 0 || range < -dmin))
	{
		valmax <- 0
		valmin <- -range
	}

	seq(from=valmin, to=valmax, by=step)
}


#' R's pretty algorithm implemented in R
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @param n number of axis intervals (specify one of \code{m} or \code{n})
#' @param min.n nonnegative integer giving the \emph{minimal} number of intervals. If \code{min.n == 0}, \code{pretty(.)} may return a single value.
#' @param shrink.sml positive numeric by a which a default scale is shrunk in the case when \code{range(x)} is very small (usually 0).
#' @param high.u.bias non-negative numeric, typically \code{> 1}. The interval unit is determined as \code{\{1,2,5,10\}} times \code{b}, a power of 10. Larger \code{high.u.bias} values favor larger units.
#' @param u5.bias non-negative numeric multiplier favoring factor 5 over 2. Default and 'optimal': \code{u5.bias = .5 + 1.5*high.u.bias}.
#' @return vector of axis label locations
#' @references
#' Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) \emph{The New S Language}. Wadsworth & Brooks/Cole.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
rpretty <- function(dmin, dmax, m=6, n=floor(m)-1, min.n=n%/%3, shrink.sml = 0.75, high.u.bias=1.5, u5.bias=0.5 + 1.5*high.u.bias)
{
	ndiv <- n
	h <- high.u.bias
	h5 <- u5.bias

	dx <- dmax-dmin
	if(dx==0 && dmax==0)
	{
		cell <- 1
		i_small <- TRUE
		U <- 1
	}
	else
	{
		cell <- max(abs(dmin), abs(dmax))
		U <- 1 + ifelse(h5 >= 1.5*h+0.5, 1/(1+h), 1.5/(1+h5))
		i_small = dx < (cell * U * max(1, ndiv) * 1e-07 * 3)
	}

	if(i_small)
	{
		if(cell > 10)
		{
			cell <- 9+cell/10
		}
    	cell <- cell * shrink.sml
		if(min.n > 1) cell <- cell/min.n
	}	
	else
	{
		cell <- dx
		if(ndiv > 1) cell <- cell/ndiv
	}

	if(cell < 20 * 1e-07)
		cell <- 20 * 1e-07
	
	base <- 10^floor(log10(cell))

	unit <- base

	if((2*base)-cell < h*(cell-unit))
	{
		unit <- 2*base
		if((5*base)-cell < h5*(cell-unit))
		{
			unit <- 5*base
			if((10*base)-cell < h*(cell-unit))
				unit <- 10*base
		}
	}

	# track down lattice labelings...

	## Maybe used to correct for the epsilon here??
	ns <- floor(dmin/unit + 1e-07)
	nu <- ceiling(dmax/unit - 1e-07)

	## Extend the range out beyond the data. Does this ever happen??
	while(ns*unit > dmin+(1e-07*unit)) ns <- ns-1
	while(nu*unit < dmax-(1e-07*unit)) nu <- nu+1


	## If we don't have quite enough labels, extend the range out to make more (these labels are beyond the data :( )
	k <- floor(0.5 + nu-ns)
	if(k < min.n)
	{
		k <- min.n - k
		if(ns >=0)
		{
			nu <- nu + k/2
			ns <- ns - k/2 + k%%2
		}
		else
		{
			ns <- ns - k/2
			nu <- nu + k/2 + k%%2
		}
		ndiv <- min.n
	}
	else
	{
		ndiv <- k
	}

	graphmin <- ns*unit
	graphmax <- nu*unit

	seq(from=graphmin, to=graphmax, by=unit)
}

#' Matplotlib's labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @return vector of axis label locations
#' @references
#' \url{http://matplotlib.sourceforge.net/}
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
matplotlib <- function(dmin, dmax, m)
{
	steps <- c(1,2,5,10)
	nbins <- m
	trim <- TRUE

	vmin <- dmin
	vmax <- dmax
	params <- .matplotlib.scale.range(vmin, vmax, nbins)
	scale <- params[1]
	offset <- params[2]

	vmin <- vmin-offset
	vmax <- vmax-offset

	rawStep <- (vmax-vmin)/nbins
	scaledRawStep <- rawStep/scale

	bestMax <- vmax
	bestMin <- vmin

	scaledStep <- 1
	chosenFactor <- 1

	for (step in steps)
	{
		if (step >= scaledRawStep)
		{
			scaledStep <- step*scale
			chosenFactor <- step
			bestMin <- scaledStep * floor(vmin/scaledStep)
			bestMax <- bestMin + scaledStep*nbins
			if (bestMax >= vmax)
				break
		}
	}
	if (trim)
	{
		extraBins <- floor((bestMax-vmax)/scaledStep)
		nbins <- nbins-extraBins
	}
	graphMin <- bestMin+offset
	graphMax <- graphMin+nbins*scaledStep

	seq(from=graphMin, to=graphMax, by=scaledStep)
}

.matplotlib.scale.range <- function(min, max, bins)
{
	threshold <- 100
	dv <- abs(max-min)
	maxabsv<-max(abs(min), abs(max))
	if (maxabsv == 0 || dv/maxabsv<10^-12)
		return(c(1, 0))

	meanv <- 0.5*(min+max)

	if ((abs(meanv)/dv) < threshold)
		offset<- 0
	else if (meanv>0)
	{
		exp<-floor(log10(meanv))
		offset = 10.0^exp
	} else
	{
		exp <- floor(log10(-1*meanv))
		offset <- -10.0^exp
	}
	exp <- floor(log10(dv/bins))
	scale = 10.0^exp
	c(scale, offset)
}



#' gnuplot's labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @return vector of axis label locations
#' @references
#' \url{http://www.gnuplot.info/}
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
gnuplot <- function(dmin, dmax, m)
{
	ntick <- floor(m)
	power <- 10^floor(log10(dmax-dmin))
	norm_range <- (dmax-dmin)/power
	p <- (ntick-1) / norm_range

	if(p > 40)
		t <- 0.05
	else if(p > 20)
		t <- 0.1
	else if(p > 10)
		t <- 0.2
	else if(p > 4)
		t <- 0.5
	else if(p > 2)
		t <- 1
	else if(p > 0.5)
		t <- 2
	else
		t <- ceiling(norm_range)

	d <- t*power
	graphmin <- floor(dmin/d) * d
	graphmax <- ceiling(dmax/d) * d

	seq(from=graphmin, to=graphmax, by=d)
}



#' Sparks' labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @return vector of axis label locations
#' @references
#' Sparks, D. N. (1971) AS 44. Scatter Diagram Plotting, Journal of the Royal Statistical Society. Series C., pp. 327-331.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
sparks <- function(dmin, dmax, m)
{
	fm <- m-1
	ratio <- 0
	key <- 1
	kount <- 0
	r <- dmax-dmin
	b <- dmin
	
	while(ratio <= 0.8)
	{
		while(key <= 2)
		{
			while(r <= 1)
			{
				kount <- kount + 1
				r <- r*10
			}
			while(r > 10)
			{
				kount <- kount - 1
				r <- r/10
			}

			b <- b*(10^kount)
			if( b < 0 && b != trunc(b)) b <- b-1
			b <- trunc(b)/(10^kount)
			r <- (dmax-b)/fm
			kount <- 0
			key <- key+2
		}
	
		fstep <- trunc(r)
		if(fstep != r) fstep <- fstep+1
		if(r < 1.5) fstep <- fstep-0.5
		fstep <- fstep/(10^kount)
		ratio <- (dmax - dmin)*(fm*fstep)
		kount <- 1
		key <- 2
	}
	fmin <- b
	c <- fstep*trunc(b/fstep)
	if(c < 0 && c != b) c <- c-fstep
	if((c+fm*fstep) > dmax) fmin <- c
	
	seq(from=fmin, to=fstep*(m-1), by=fstep)
}


#' Thayer and Storer's labeling algorithm
#'
#' @param dmin minimum of the data range
#' @param dmax maximum of the data range
#' @param m number of axis labels
#' @return vector of axis label locations
#' @references
#' Thayer, R. P. and Storer, R. F. (1969) AS 21. Scale Selection for Computer Plots, Journal of the Royal Statistical Society. Series C., pp. 206-208.
#' @author Justin Talbot \email{jtalbot@@stanford.edu}
#' @export
thayer <- function(dmin, dmax, m)
{
	r <- dmax-dmin
	b <- dmin
	kount <- 0
	kod <- 0

	while(kod < 2)
	{
		while(r <= 1)
		{
			kount <- kount+1
			r <- r*10
		}
		while(r > 10)
		{
			kount <- kount-1
			r <- r/10
		}
		b <- b*(10^kount)
		if(b < 0)
			b <- b-1
		ib <- trunc(b)
		b <- ib
		b <- b/(10^kount)
		r <- dmax-b
		a <- r/(m-1)
		kount <- 0
		while(a <= 1)
		{
			kount <- kount+1
			a <- a*10
		}
		while(a > 10)
		{
			kount <- kount-1
			a <- a/10
		}
		ia <- trunc(a)
		if(ia == 6) ia <- 7
		if(ia == 8) ia <- 9

		aa <- 0
		if(a < 1.5) aa <- -0.5
		a <- aa + 1 + ia
		a <- a/(10^kount)
				
		test <- (m-1) * a
		test1 <- (dmax-dmin)/test
		if(test1 > 0.8)
			kod <- 2

		if(kod < 2)
		{
			kount <- 1
			r <- dmax-dmin
			b <- dmin
			kod <- kod + 1
		}
	}

	iab <- trunc(b/a)
	if(iab < 0) iab <- iab-1
	c <- a * iab
	d <- c + (m-1)*a
	if(d >= dmax)
		b <- c

	valmin <- b
	valmax <- b + a*(m-1)	

	seq(from=valmin, to=valmax, by=a)
}