#######################################################
#                                                     #
# S-plus functions for Minnotte's Local Mode Test     #
#                                                     #
# Save as text, then source into S-plus               #
#                                                     #
# Copyright 1995, Michael C. Minnotte                 #
#                                                     #
# Requires fortran code for calculations              #
#                                                     #
# For assistance, or to report bugs, please contact:  #
# minnotte@math.usu.edu                               #
#                                                     #
# For more information on the test, see               #
# Minnotte, M.C., (1997), "Nonparametric testing of   #
# the existence of modes," The Annals of Statistics,  #
# 25, 1646-1660.                                      #
#                                                     #
# Functions include:                                  #
#                                                     #
# tme.f - provides input and output for fortran       #
#    code.  Assumes tme executable is available       #
#    in directory S-plus is called from.              #
# treeplot - plots a mode tree, including indicators  #
#    of mode test significance                        #
# nicerange - calculates an extension of the range of #
#    the data for density estimation and mode tree    #
#    purposes.                                        #
#                                                     #
#######################################################



tme.f<-function(data,doplot=F,maintit="")
{
# Calls mode testing executable tme.
# Assumes tme sits in directory S-plus was called from
# data - univariate vector of potentially multimodal data
# doplot - if T, produces a mode tree at the end of the run indicating
#   significance of modes at the 0.05 level
# maintit - overall title for mode tree from doplot
	write(data,"dataxyz",1,F)
	unix("rm trackingxyz")
	unix("tme > trackingxyz")
	tree<-matrix(scan("treeoutxyz"),ncol=4,byrow=T)
	modes<-matrix(scan("modeoutxyz"),ncol=5,byrow=T)
#	unix("rm dataxyz treeoutxyz modeoutxyz")
# uncomment to clean up directory after run
        out <- list(tree = tree, modes = modes)
        if(doplot)
                treeplot(out, maintit)
        return(out)
}

treeplot<-
function(modetree, maintit = "", sp = T, sl = 0.05, nm = 100)
{
# mode tree plotting function
# modetree - list consisting of two matrices tree and modes, as from
#   tme.f
# sp - show significant p-values, if T, plot full octagons for modes
#   tested with p < sl, empty octagons for those with p >= sl
# sl - significance level break point
# nm - maximum number of modes to show
	tree <- modetree$tree
	modes <- modetree$modes
	if(nm < nrow(modes)) {
		modes <- modes[1:nm,  ]
		tree <- tree[tree[, 2] >= modes[nm, 2],  ]
	}
	nr <- nrow(tree)
	plot(tree[floor(nr/2), 1], tree[floor(nr/2), 2], type = "n", xlim = 
		nicerange(tree[, 1]), ylim = range(tree[, 2]), log =
		"y", xlab = "",ylab = "h", main = maintit)
	t1 <- (1:nrow(tree))[tree[, 4] == 1]
	t2 <- (1:nrow(tree))[tree[, 4] == 2]
	segments(tree[t1, 1], tree[t1, 2], tree[tree[t1, 3], 1], tree[tree[
		t1, 3], 2], lty = 1)
	segments(tree[t2, 1], tree[t2, 2], tree[tree[t2, 3], 1], tree[tree[
		t2, 3], 2], lty = 2)
	if(sp) {
		good <- (1:nrow(modes))[modes[, 3] < sl]
		if(length(good) > 0)
			points(modes[good, 1], modes[good, 2], pch = 16)
		if(length(good) < nrow(modes))
			points(modes[ - good, 1], modes[ - good, 2], pch = 1)
	}
}


nicerange<-function(x, beta = 0.2)
{
#  extends range of data (x) by beta/2*range on each end
        ab <- range(x)  # interval surrounding data
        del <- ((ab[2] - ab[1]) * beta)/2
        return(c(ab + c( - del, del)))
}